MOAB: Mesh Oriented datABase  (version 5.2.1)
FindPtFuncs.h
Go to the documentation of this file.
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines