MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 #ifndef FINDPTFUNCS_H 00002 #define FINDPTFUNCS_H 00003 00004 #include "float.h" 00005 #include "math.h" 00006 00007 /*====================================================== 00008 / from types.h 00009 /====================================================== 00010 */ 00011 00012 #define MOAB_POLY_EPS (128*DBL_EPSILON) 00013 #define MOAB_POLY_PI 3.1415926535897932384626433832795028841971693993751058209749445923 00014 00015 #define mbsqrt sqrt 00016 #define mbabs fabs 00017 #define mbcos cos 00018 #define mbsin sin 00019 #define mbfloor floor 00020 #define mbceil ceil 00021 00022 #ifdef INTEGER 00023 #undef INTEGER 00024 #endif 00025 00026 #ifdef real 00027 #undef real 00028 #endif 00029 00030 /* integer type to use for everything */ 00031 #if defined(USE_LONG) 00032 # define INTEGER long 00033 #elif defined(USE_LONG_LONG) 00034 # define INTEGER long long 00035 #elif defined(USE_SHORT) 00036 # define INTEGER short 00037 #else 00038 # define INTEGER int 00039 #endif 00040 00041 /* when defined, use the given type for global indices instead of INTEGER */ 00042 #if defined(USE_GLOBAL_LONG_LONG) 00043 # define GLOBAL_INT long long 00044 #elif defined(USE_GLOBAL_LONG) 00045 # define GLOBAL_INT long 00046 #else 00047 # define GLOBAL_INT long 00048 #endif 00049 00050 /* floating point type to use for everything */ 00051 #if defined(USE_FLOAT) 00052 typedef float realType; 00053 #elif defined(USE_LONG_DOUBLE) 00054 typedef long double realType; 00055 #else 00056 typedef double realType; 00057 #endif 00058 00059 /* apparently uint and ulong can be defined already in standard headers */ 00060 #ifndef uint 00061 typedef signed INTEGER sint; 00062 #endif 00063 00064 #ifndef uint 00065 typedef unsigned INTEGER uint; 00066 #endif 00067 #undef INTEGER 00068 00069 #ifndef slong 00070 #ifdef GLOBAL_INT 00071 typedef signed GLOBAL_INT slong; 00072 #else 00073 typedef sint slong; 00074 #endif 00075 #endif 00076 00077 #ifndef ulong 00078 #ifdef GLOBAL_INT 00079 typedef unsigned GLOBAL_INT ulong; 00080 #else 00081 typedef uint ulong; 00082 #endif 00083 #endif 00084 00085 /*====================================================== 00086 / from poly.h 00087 /====================================================== 00088 */ 00089 00090 /* 00091 For brevity's sake, some names have been shortened 00092 Quadrature rules 00093 Gauss -> Gauss-Legendre quadrature (open) 00094 Lobatto -> Gauss-Lobatto-Legendre quadrature (closed at both ends) 00095 Polynomial bases 00096 Legendre -> Legendre basis 00097 Gauss -> Lagrangian basis using Gauss quadrature nodes 00098 Lobatto -> Lagrangian basis using Lobatto quadrature nodes 00099 */ 00100 00101 /*-------------------------------------------------------------------------- 00102 Legendre Polynomial Matrix Computation 00103 (compute P_i(x_j) for i = 0, ..., n and a given set of x) 00104 --------------------------------------------------------------------------*/ 00105 00106 /* precondition: n >= 1 00107 inner index is x index (0 ... m-1); 00108 outer index is Legendre polynomial number (0 ... n) 00109 */ 00110 void legendre_matrix(const realType *x, int m, realType *P, int n); 00111 00112 /* precondition: n >= 1 00113 inner index is Legendre polynomial number (0 ... n) 00114 outer index is x index (0 ... m-1); 00115 */ 00116 void legendre_matrix_t(const realType *x, int m, realType *P, int n); 00117 00118 /* precondition: n >= 1 00119 compute P_i(x) with i = 0 ... n 00120 */ 00121 void legendre_row(realType x, realType *P, int n); 00122 00123 00124 /*-------------------------------------------------------------------------- 00125 Quadrature Nodes and Weights Calculation 00126 00127 call the _nodes function before calling the _weights function 00128 --------------------------------------------------------------------------*/ 00129 00130 void gauss_nodes(realType *z, int n); /* n nodes (order = 2n-1) */ 00131 void lobatto_nodes(realType *z, int n); /* n nodes (order = 2n-3) */ 00132 00133 void gauss_weights(const realType *z, realType *w, int n); 00134 void lobatto_weights(const realType *z, realType *w, int n); 00135 00136 /*-------------------------------------------------------------------------- 00137 Lagrangian to Legendre Change-of-basis Matrix 00138 --------------------------------------------------------------------------*/ 00139 00140 /* precondition: n >= 2 00141 given the Gauss quadrature rule (z,w,n), compute the square matrix J 00142 for transforming from the Gauss basis to the Legendre basis: 00143 00144 u_legendre(i) = sum_j J(i,j) u_gauss(j) 00145 00146 computes J = .5 (2i+1) w P (z ) 00147 ij j i j 00148 00149 in column major format (inner index is i, the Legendre index) 00150 */ 00151 void gauss_to_legendre(const realType *z, const realType *w, int n, realType *J); 00152 00153 /* precondition: n >= 2 00154 same as above, but 00155 in row major format (inner index is j, the Gauss index) 00156 */ 00157 void gauss_to_legendre_t(const realType *z, const realType *w, int n, realType *J); 00158 00159 /* precondition: n >= 3 00160 given the Lobatto quadrature rule (z,w,n), compute the square matrix J 00161 for transforming from the Lobatto basis to the Legendre basis: 00162 00163 u_legendre(i) = sum_j J(i,j) u_lobatto(j) 00164 00165 in column major format (inner index is i, the Legendre index) 00166 */ 00167 void lobatto_to_legendre(const realType *z, const realType *w, int n, realType *J); 00168 00169 /*-------------------------------------------------------------------------- 00170 Lagrangian basis function evaluation 00171 --------------------------------------------------------------------------*/ 00172 00173 /* given the Lagrangian nodes (z,n) and evaluation points (x,m) 00174 evaluate all Lagrangian basis functions at all points x 00175 00176 inner index of output J is the basis function index (row-major format) 00177 provide work array with space for 4*n doubles 00178 */ 00179 void lagrange_weights(const realType *z, unsigned n, 00180 const realType *x, unsigned m, 00181 realType *J, realType *work); 00182 00183 /* given the Lagrangian nodes (z,n) and evaluation points (x,m) 00184 evaluate all Lagrangian basis functions and their derivatives 00185 00186 inner index of outputs J,D is the basis function index (row-major format) 00187 provide work array with space for 6*n doubles 00188 */ 00189 void lagrange_weights_deriv(const realType *z, unsigned n, 00190 const realType *x, unsigned m, 00191 realType *J, realType *D, realType *work); 00192 00193 /*-------------------------------------------------------------------------- 00194 Speedy Lagrangian Interpolation 00195 00196 Usage: 00197 00198 lagrange_data p; 00199 lagrange_setup(&p,z,n); * setup for nodes z[0 ... n-1] * 00200 00201 the weights 00202 p->J [0 ... n-1] interpolation weights 00203 p->D [0 ... n-1] 1st derivative weights 00204 p->D2[0 ... n-1] 2nd derivative weights 00205 are computed for a given x with: 00206 lagrange_0(p,x); * compute p->J * 00207 lagrange_1(p,x); * compute p->J, p->D * 00208 lagrange_2(p,x); * compute p->J, p->D, p->D2 * 00209 lagrange_2u(p); * compute p->D2 after call of lagrange_1(p,x); * 00210 These functions use the z array supplied to setup 00211 (that pointer should not be freed between calls) 00212 Weights for x=z[0] and x=z[n-1] are computed during setup; access as: 00213 p->J_z0, etc. and p->J_zn, etc. 00214 00215 lagrange_free(&p); * deallocate memory allocated by setup * 00216 --------------------------------------------------------------------------*/ 00217 00218 typedef struct { 00219 unsigned n; /* number of Lagrange nodes */ 00220 const realType *z; /* Lagrange nodes (user-supplied) */ 00221 realType *J, *D, *D2; /* weights for 0th,1st,2nd derivatives */ 00222 realType *J_z0, *D_z0, *D2_z0; /* ditto at z[0] (computed at setup) */ 00223 realType *J_zn, *D_zn, *D2_zn; /* ditto at z[n-1] (computed at setup) */ 00224 realType *w, *d, *u0, *v0, *u1, *v1, *u2, *v2; /* work data */ 00225 } lagrange_data; 00226 00227 void lagrange_setup(lagrange_data *p, const realType *z, unsigned n); 00228 void lagrange_free(lagrange_data *p); 00229 00230 void lagrange_0(lagrange_data *p, realType x) ; 00231 void lagrange_1(lagrange_data *p, realType x) ; 00232 void lagrange_2(lagrange_data *p, realType x) ; 00233 void lagrange_2u(lagrange_data *p) ; 00234 00235 /*====================================================== 00236 / from tensor.h 00237 /====================================================== 00238 */ 00239 00240 /*-------------------------------------------------------------------------- 00241 1-,2-,3-d Tensor Application 00242 00243 the 3d case: 00244 tensor_f3(R,mr,nr, S,ms,ns, T,mt,nt, u,v, work1,work2) 00245 gives v = [ R (x) S (x) T ] u 00246 where R is mr x nr, S is ms x ns, T is mt x nt, 00247 each in row- or column-major format according to f := r | c 00248 u is nr x ns x nt in column-major format (inner index is r) 00249 v is mr x ms x mt in column-major format (inner index is r) 00250 --------------------------------------------------------------------------*/ 00251 00252 void tensor_c1(const realType *R, unsigned mr, unsigned nr, 00253 const realType *u, realType *v); 00254 void tensor_r1(const realType *R, unsigned mr, unsigned nr, 00255 const realType *u, realType *v); 00256 00257 /* work holds mr*ns realTypes */ 00258 void tensor_c2(const realType *R, unsigned mr, unsigned nr, 00259 const realType *S, unsigned ms, unsigned ns, 00260 const realType *u, realType *v, realType *work); 00261 void tensor_r2(const realType *R, unsigned mr, unsigned nr, 00262 const realType *S, unsigned ms, unsigned ns, 00263 const realType *u, realType *v, realType *work); 00264 00265 /* work1 holds mr*ns*nt realTypes, 00266 work2 holds mr*ms*nt realTypes */ 00267 void tensor_c3(const realType *R, unsigned mr, unsigned nr, 00268 const realType *S, unsigned ms, unsigned ns, 00269 const realType *T, unsigned mt, unsigned nt, 00270 const realType *u, realType *v, realType *work1, realType *work2); 00271 void tensor_r3(const realType *R, unsigned mr, unsigned nr, 00272 const realType *S, unsigned ms, unsigned ns, 00273 const realType *T, unsigned mt, unsigned nt, 00274 const realType *u, realType *v, realType *work1, realType *work2); 00275 00276 /*-------------------------------------------------------------------------- 00277 1-,2-,3-d Tensor Application of Row Vectors (for Interpolation) 00278 00279 the 3d case: 00280 v = tensor_i3(Jr,nr, Js,ns, Jt,nt, u, work) 00281 same effect as tensor_r3(Jr,1,nr, Js,1,ns, Jt,1,nt, u,&v, work1,work2): 00282 gives v = [ Jr (x) Js (x) Jt ] u 00283 where Jr, Js, Jt are row vectors (interpolation weights) 00284 u is nr x ns x nt in column-major format (inner index is r) 00285 v is a scalar 00286 --------------------------------------------------------------------------*/ 00287 00288 realType tensor_i1(const realType *Jr, unsigned nr, const realType *u); 00289 00290 /* work holds ns realTypes */ 00291 realType tensor_i2(const realType *Jr, unsigned nr, 00292 const realType *Js, unsigned ns, 00293 const realType *u, realType *work); 00294 00295 /* work holds ns*nt + nt realTypes */ 00296 realType tensor_i3(const realType *Jr, unsigned nr, 00297 const realType *Js, unsigned ns, 00298 const realType *Jt, unsigned nt, 00299 const realType *u, realType *work); 00300 00301 /*-------------------------------------------------------------------------- 00302 1-,2-,3-d Tensor Application of Row Vectors 00303 for simultaneous Interpolation and Gradient computation 00304 00305 the 3d case: 00306 v = tensor_ig3(Jr,Dr,nr, Js,Ds,ns, Jt,Dt,nt, u,g, work) 00307 gives v = [ Jr (x) Js (x) Jt ] u 00308 g_0 = [ Dr (x) Js (x) Jt ] u 00309 g_1 = [ Jr (x) Ds (x) Jt ] u 00310 g_2 = [ Jr (x) Js (x) Dt ] u 00311 where Jr,Dr,Js,Ds,Jt,Dt are row vectors 00312 (interpolation & derivative weights) 00313 u is nr x ns x nt in column-major format (inner index is r) 00314 v is a scalar, g is an array of 3 realTypes 00315 --------------------------------------------------------------------------*/ 00316 00317 realType tensor_ig1(const realType *Jr, const realType *Dr, unsigned nr, 00318 const realType *u, realType *g); 00319 00320 /* work holds 2*ns realTypes */ 00321 realType tensor_ig2(const realType *Jr, const realType *Dr, unsigned nr, 00322 const realType *Js, const realType *Ds, unsigned ns, 00323 const realType *u, realType *g, realType *work); 00324 00325 /* work holds 2*ns*nt + 3*ns realTypes */ 00326 realType tensor_ig3(const realType *Jr, const realType *Dr, unsigned nr, 00327 const realType *Js, const realType *Ds, unsigned ns, 00328 const realType *Jt, const realType *Dt, unsigned nt, 00329 const realType *u, realType *g, realType *work); 00330 00331 /*====================================================== 00332 / from findpt.h 00333 /====================================================== 00334 */ 00335 00336 typedef struct { 00337 realType x[2], A[4], axis_bnd[4]; 00338 } obbox_2; 00339 00340 typedef struct { 00341 realType x[3], A[9], axis_bnd[6]; 00342 } obbox_3; 00343 00344 typedef struct { 00345 unsigned hash_n; 00346 realType bnd[4]; /* bounds for all elements */ 00347 realType fac[2]; /* fac[i] = hash_n / (bnd[2*i+1]-bnd[2*i]) */ 00348 obbox_2 *obb; /* obb[nel] -- bounding box info for each element */ 00349 uint *offset; /* hash table -- for cell i,j: 00350 uint index = j*hash_n+i, 00351 b = offset[index ], 00352 e = offset[index+1]; 00353 elements in cell are 00354 offset[b], offset[b+1], ..., offset[e-1] */ 00355 unsigned max; /* maximum # of elements in any cell */ 00356 } hash_data_2; 00357 00358 typedef struct { 00359 unsigned hash_n; 00360 realType bnd[6]; /* bounds for all elements */ 00361 realType fac[3]; /* fac[i] = hash_n / (bnd[2*i+1]-bnd[2*i]) */ 00362 obbox_3 *obb; /* obb[nel] -- bounding box info for each element */ 00363 uint *offset; /* hash table -- for cell i,j,k: 00364 uint index = (k*hash_n+j)*hash_n+i, 00365 b = offset[index ], 00366 e = offset[index+1]; 00367 elements in cell are 00368 offset[b], offset[b+1], ..., offset[e-1] */ 00369 unsigned max; /* maximum # of elements in any cell */ 00370 } hash_data_3; 00371 00372 typedef struct { 00373 uint el; 00374 realType r[3]; 00375 realType dist; 00376 } findpt_listel; 00377 00378 typedef struct { 00379 unsigned constraints; 00380 unsigned de, d1; 00381 realType *x[2], *fd1[2]; 00382 } opt_edge_data_2; 00383 00384 typedef struct { 00385 unsigned constraints; 00386 realType x[2], jac[4]; 00387 } opt_point_data_2; 00388 00389 typedef struct { 00390 lagrange_data *ld; 00391 unsigned size[3]; 00392 const realType *elx[2]; 00393 opt_edge_data_2 ed; 00394 opt_point_data_2 pd; 00395 realType *work; 00396 realType x[2], jac[4]; 00397 } opt_data_2; 00398 00399 typedef struct { 00400 const realType *xw[2]; /* geometry data */ 00401 realType *z[2]; /* lobatto nodes */ 00402 lagrange_data ld[2]; /* interpolation, derivative weights & data */ 00403 unsigned nptel; /* nodes per element */ 00404 hash_data_2 *hash; /* geometric hashing data */ 00405 findpt_listel *list, **sorted, **end; 00406 /* pre-allocated list of elements to 00407 check (found by hashing), and 00408 pre-allocated list of pointers into 00409 the first list for sorting */ 00410 opt_data_2 *od; /* data for the optimization algorithm */ 00411 realType *od_work; 00412 } findpt_data_2; 00413 00414 typedef struct { 00415 unsigned constraints; 00416 unsigned dn, d1, d2; 00417 realType *x[3], *fdn[3]; 00418 } opt_face_data_3; 00419 00420 typedef struct { 00421 unsigned constraints; 00422 unsigned de, d1, d2; 00423 realType *x[3], *fd1[3], *fd2[3]; 00424 } opt_edge_data_3; 00425 00426 typedef struct { 00427 unsigned constraints; 00428 realType x[3], jac[9]; 00429 } opt_point_data_3; 00430 00431 typedef struct { 00432 lagrange_data *ld; 00433 unsigned size[4]; 00434 const realType *elx[3]; 00435 opt_face_data_3 fd; 00436 opt_edge_data_3 ed; 00437 opt_point_data_3 pd; 00438 realType *work; 00439 realType x[3], jac[9]; 00440 } opt_data_3; 00441 00442 typedef struct { 00443 const realType *xw[3]; /* geometry data */ 00444 realType *z[3]; /* lobatto nodes */ 00445 lagrange_data ld[3]; /* interpolation, derivative weights & data */ 00446 unsigned nptel; /* nodes per element */ 00447 hash_data_3 *hash; /* geometric hashing data */ 00448 findpt_listel *list, **sorted, **end; 00449 /* pre-allocated list of elements to 00450 check (found by hashing), and 00451 pre-allocated list of pointers into 00452 the first list for sorting */ 00453 opt_data_3 *od; /* data for the optimization algorithm */ 00454 realType *od_work; 00455 } findpt_data_3; 00456 00457 findpt_data_2 *findpt_setup_2( 00458 const realType *const xw[2], const unsigned n[2], uint nel, 00459 uint max_hash_size, realType bbox_tol); 00460 findpt_data_3 *findpt_setup_3( 00461 const realType *const xw[3], const unsigned n[3], uint nel, 00462 uint max_hash_size, realType bbox_tol); 00463 00464 void findpt_free_2(findpt_data_2 *p); 00465 void findpt_free_3(findpt_data_3 *p); 00466 00467 const realType *findpt_allbnd_2(const findpt_data_2 *p); 00468 const realType *findpt_allbnd_3(const findpt_data_3 *p); 00469 00470 typedef int (*findpt_func)(void *, const realType *, int, uint *, realType *, realType *); 00471 int findpt_2(findpt_data_2 *p, const realType x[2], int guess, 00472 uint *el, realType r[2], realType *dist); 00473 int findpt_3(findpt_data_3 *p, const realType x[3], int guess, 00474 uint *el, realType r[3], realType *dist); 00475 00476 void findpt_weights_2(findpt_data_2 *p, const realType r[2]); 00477 00478 void findpt_weights_3(findpt_data_3 *p, const realType r[3]); 00479 00480 double findpt_eval_2(findpt_data_2 *p, const realType *u); 00481 00482 double findpt_eval_3(findpt_data_3 *p, const realType *u); 00483 00484 /*====================================================== 00485 / from extrafindpt.h 00486 /====================================================== 00487 */ 00488 00489 void opt_alloc_3(opt_data_3 *p, lagrange_data *ld); 00490 void opt_free_3(opt_data_3 *p); 00491 double opt_findpt_3(opt_data_3 *p, const realType *const elx[3], 00492 const realType xstar[3], realType r[3], unsigned *constr); 00493 void opt_vol_set_intp_3(opt_data_3 *p, const realType r[3]); 00494 00495 /* for 2d spectralQuad */ 00496 /*-------------------------------------------------------------------------- 00497 00498 2 - D 00499 00500 --------------------------------------------------------------------------*/ 00501 00502 void opt_alloc_2(opt_data_2 *p, lagrange_data *ld); 00503 void opt_free_2(opt_data_2 *p); 00504 double opt_findpt_2(opt_data_2 *p, const realType *const elx[2], 00505 const realType xstar[2], realType r[2], unsigned *constr); 00506 00507 extern const unsigned opt_no_constraints_2; 00508 extern const unsigned opt_no_constraints_3; 00509 00510 /*====================================================== 00511 / from errmem.h 00512 /====================================================== 00513 */ 00514 00515 /* requires: 00516 <stdlib.h> for malloc, calloc, realloc, free 00517 */ 00518 00519 /*-------------------------------------------------------------------------- 00520 Error Reporting 00521 Memory Allocation Wrappers to Catch Out-of-memory 00522 --------------------------------------------------------------------------*/ 00523 00524 /* #include "malloc.h" */ 00525 #include <stdlib.h> 00526 00527 #ifdef __GNUC__ 00528 void fail(const char *fmt, ...) __attribute__ ((noreturn)); 00529 #define MAYBE_UNUSED __attribute__ ((unused)) 00530 #else 00531 void fail(const char *fmt, ...); 00532 #define MAYBE_UNUSED 00533 #endif 00534 00535 #if 0 00536 {} 00537 #endif 00538 00539 static void *smalloc(size_t size, const char *file) MAYBE_UNUSED; 00540 static void *smalloc(size_t size, const char *file) 00541 { 00542 void *res = malloc(size); 00543 if(!res && size) fail("%s: allocation of %d bytes failed\n",file,(int)size); 00544 return res; 00545 } 00546 00547 static void *scalloc(size_t nmemb, size_t size, const char *file) MAYBE_UNUSED; 00548 static void *scalloc(size_t nmemb, size_t size, const char *file) 00549 { 00550 void *res = calloc(nmemb, size); 00551 if(!res && nmemb) 00552 fail("%s: allocation of %d bytes failed\n",file,(int)size*nmemb); 00553 return res; 00554 } 00555 00556 static void *srealloc(void *ptr, size_t size, const char *file) MAYBE_UNUSED; 00557 static void *srealloc(void *ptr, size_t size, const char *file) 00558 { 00559 void *res = realloc(ptr, size); 00560 if(!res && size) fail("%s: allocation of %d bytes failed\n",file,(int)size); 00561 return res; 00562 } 00563 00564 #define tmalloc(type, count) \ 00565 ((type*) smalloc((count)*sizeof(type),__FILE__) ) 00566 #define tcalloc(type, count) \ 00567 ((type*) scalloc((count),sizeof(type),__FILE__) ) 00568 #define trealloc(type, ptr, count) \ 00569 ((type*) srealloc((ptr),(count)*sizeof(type),__FILE__) ) 00570 /* 00571 typedef struct { size_t size; void *ptr; } buffer; 00572 00573 static void buffer_init_(buffer *b, size_t size, const char *file) MAYBE_UNUSED; 00574 static void buffer_init_(buffer *b, size_t size, const char *file) 00575 { 00576 b->size=size, b->ptr=smalloc(size,file); 00577 } 00578 static void buffer_reserve_(buffer *b, size_t min, const char *file) 00579 MAYBE_UNUSED; 00580 static void buffer_reserve_(buffer *b, size_t min, const char *file) 00581 { 00582 size_t size = b->size; 00583 if(size<min) { 00584 size+=size/2+1; 00585 if(size<min) size=min; 00586 b->ptr=srealloc(b->ptr,size,file); 00587 } 00588 } 00589 static void buffer_free(buffer *b) MAYBE_UNUSED; 00590 static void buffer_free(buffer *b) { free(b->ptr); } 00591 00592 #define buffer_init(b,size) buffer_init_(b,size,__FILE__) 00593 #define buffer_reserve(b,min) buffer_reserve_(b,min,__FILE__) 00594 */ 00595 00596 00597 /*====================================================== 00598 / from minmax.h 00599 /====================================================== 00600 */ 00601 00602 /*-------------------------------------------------------------------------- 00603 Min, Max, Norm 00604 --------------------------------------------------------------------------*/ 00605 00606 #ifdef __GNUC__ 00607 #define MAYBE_UNUSED __attribute__ ((unused)) 00608 #else 00609 #define MAYBE_UNUSED 00610 #endif 00611 00612 #define DECLMINMAX(type, prefix) \ 00613 static type prefix##min_2(type a, type b) MAYBE_UNUSED; \ 00614 static type prefix##min_2(type a, type b) { \ 00615 return b<a?b:a; \ 00616 } \ 00617 static type prefix##max_2(type a, type b) MAYBE_UNUSED; \ 00618 static type prefix##max_2(type a, type b) { \ 00619 return a>b?a:b; \ 00620 } \ 00621 static void prefix##minmax_2(type *min, type *max, type a, \ 00622 type b) MAYBE_UNUSED; \ 00623 static void prefix##minmax_2(type *min, type *max, type a, type b) { \ 00624 if(b<a) *min=b, *max=a; \ 00625 else *min=a, *max=b; \ 00626 } \ 00627 static type prefix##min_3(type a, type b, type c) MAYBE_UNUSED; \ 00628 static type prefix##min_3(type a, type b, type c) { \ 00629 return b<a?(c<b?c:b):(c<a?c:a); \ 00630 } \ 00631 static type prefix##max_3(type a, type b, type c) MAYBE_UNUSED; \ 00632 static type prefix##max_3(type a, type b, type c) { \ 00633 return a>b?(a>c?a:c):(b>c?b:c); \ 00634 } \ 00635 static void prefix##minmax_3(type *min, type *max, type a, type b, \ 00636 type c) MAYBE_UNUSED; \ 00637 static void prefix##minmax_3(type *min, type *max, type a, type b, \ 00638 type c) { \ 00639 if(b<a) *min=prefix##min_2(b,c), *max=prefix##max_2(a,c); \ 00640 else *min=prefix##min_2(a,c), *max=prefix##max_2(b,c); \ 00641 } 00642 00643 DECLMINMAX(int, i) 00644 DECLMINMAX(unsigned, u) 00645 DECLMINMAX(realType, r) 00646 #undef DECLMINMAX 00647 00648 static realType r1norm_1(realType a) MAYBE_UNUSED; 00649 static realType r1norm_1(realType a) { 00650 return mbabs(a); 00651 } 00652 static realType r1norm_2(realType a, realType b) MAYBE_UNUSED; 00653 static realType r1norm_2(realType a, realType b) { 00654 return mbabs(a)+mbabs(b); 00655 } 00656 static realType r1norm_3(realType a, realType b, realType c) MAYBE_UNUSED; 00657 static realType r1norm_3(realType a, realType b, realType c) { 00658 return mbabs(a)+mbabs(b)+mbabs(c); 00659 } 00660 static realType r2norm_1(realType a) MAYBE_UNUSED; 00661 static realType r2norm_1(realType a) { 00662 return mbsqrt(a*a); 00663 } 00664 static realType r2norm_2(realType a, realType b) MAYBE_UNUSED; 00665 static realType r2norm_2(realType a, realType b) { 00666 return mbsqrt(a*a+b*b); 00667 } 00668 static realType r2norm_3(realType a, realType b, realType c) MAYBE_UNUSED; 00669 static realType r2norm_3(realType a, realType b, realType c) { 00670 return mbsqrt(a*a+b*b+c*c); 00671 } 00672 00673 #endif 00674 00675