Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
FindPtFuncs.h
Go to the documentation of this file.
1 #ifndef FINDPTFUNCS_H
2 #define FINDPTFUNCS_H
3 
4 #include "float.h"
5 #include "math.h"
6 
7 /*======================================================
8 / from types.h
9 /======================================================
10 */
11 
12 #define MOAB_POLY_EPS ( 128 * DBL_EPSILON )
13 #define MOAB_POLY_PI 3.1415926535897932384626433832795028841971693993751058209749445923
14 
15 #define mbsqrt sqrt
16 #define mbabs fabs
17 #define mbcos cos
18 #define mbsin sin
19 #define mbfloor floor
20 #define mbceil ceil
21 
22 #ifdef INTEGER
23 #undef INTEGER
24 #endif
25 
26 #ifdef real
27 #undef real
28 #endif
29 
30 /* integer type to use for everything */
31 #if defined( USE_LONG )
32 #define INTEGER long
33 #elif defined( USE_LONG_LONG )
34 #define INTEGER long long
35 #elif defined( USE_SHORT )
36 #define INTEGER short
37 #else
38 #define INTEGER int
39 #endif
40 
41 /* when defined, use the given type for global indices instead of INTEGER */
42 #if defined( USE_GLOBAL_LONG_LONG )
43 #define GLOBAL_INT long long
44 #elif defined( USE_GLOBAL_LONG )
45 #define GLOBAL_INT long
46 #else
47 #define GLOBAL_INT long
48 #endif
49 
50 /* floating point type to use for everything */
51 #if defined( USE_FLOAT )
52 typedef float realType;
53 #elif defined( USE_LONG_DOUBLE )
54 typedef long double realType;
55 #else
56 typedef double realType;
57 #endif
58 
59 /* apparently uint and ulong can be defined already in standard headers */
60 #ifndef uint
61 typedef signed INTEGER sint;
62 #endif
63 
64 #ifndef uint
65 typedef unsigned INTEGER uint;
66 #endif
67 #undef INTEGER
68 
69 #ifndef slong
70 #ifdef GLOBAL_INT
71 typedef signed GLOBAL_INT slong;
72 #else
73 typedef sint slong;
74 #endif
75 #endif
76 
77 #ifndef ulong
78 #ifdef GLOBAL_INT
79 typedef unsigned GLOBAL_INT ulong;
80 #else
81 typedef uint ulong;
82 #endif
83 #endif
84 
85 /*======================================================
86 / from poly.h
87 /======================================================
88 */
89 
90 /*
91  For brevity's sake, some names have been shortened
92  Quadrature rules
93  Gauss -> Gauss-Legendre quadrature (open)
94  Lobatto -> Gauss-Lobatto-Legendre quadrature (closed at both ends)
95  Polynomial bases
96  Legendre -> Legendre basis
97  Gauss -> Lagrangian basis using Gauss quadrature nodes
98  Lobatto -> Lagrangian basis using Lobatto quadrature nodes
99 */
100 
101 /*--------------------------------------------------------------------------
102  Legendre Polynomial Matrix Computation
103  (compute P_i(x_j) for i = 0, ..., n and a given set of x)
104  --------------------------------------------------------------------------*/
105 
106 /* precondition: n >= 1
107  inner index is x index (0 ... m-1);
108  outer index is Legendre polynomial number (0 ... n)
109  */
110 void legendre_matrix( const realType* x, int m, realType* P, int n );
111 
112 /* precondition: n >= 1
113  inner index is Legendre polynomial number (0 ... n)
114  outer index is x index (0 ... m-1);
115  */
116 void legendre_matrix_t( const realType* x, int m, realType* P, int n );
117 
118 /* precondition: n >= 1
119  compute P_i(x) with i = 0 ... n
120  */
121 void legendre_row( realType x, realType* P, int n );
122 
123 /*--------------------------------------------------------------------------
124  Quadrature Nodes and Weights Calculation
125 
126  call the _nodes function before calling the _weights function
127  --------------------------------------------------------------------------*/
128 
129 void gauss_nodes( realType* z, int n ); /* n nodes (order = 2n-1) */
130 void lobatto_nodes( realType* z, int n ); /* n nodes (order = 2n-3) */
131 
132 void gauss_weights( const realType* z, realType* w, int n );
133 void lobatto_weights( const realType* z, realType* w, int n );
134 
135 /*--------------------------------------------------------------------------
136  Lagrangian to Legendre Change-of-basis Matrix
137  --------------------------------------------------------------------------*/
138 
139 /* precondition: n >= 2
140  given the Gauss quadrature rule (z,w,n), compute the square matrix J
141  for transforming from the Gauss basis to the Legendre basis:
142 
143  u_legendre(i) = sum_j J(i,j) u_gauss(j)
144 
145  computes J = .5 (2i+1) w P (z )
146  ij j i j
147 
148  in column major format (inner index is i, the Legendre index)
149  */
150 void gauss_to_legendre( const realType* z, const realType* w, int n, realType* J );
151 
152 /* precondition: n >= 2
153  same as above, but
154  in row major format (inner index is j, the Gauss index)
155  */
156 void gauss_to_legendre_t( const realType* z, const realType* w, int n, realType* J );
157 
158 /* precondition: n >= 3
159  given the Lobatto quadrature rule (z,w,n), compute the square matrix J
160  for transforming from the Lobatto basis to the Legendre basis:
161 
162  u_legendre(i) = sum_j J(i,j) u_lobatto(j)
163 
164  in column major format (inner index is i, the Legendre index)
165  */
166 void lobatto_to_legendre( const realType* z, const realType* w, int n, realType* J );
167 
168 /*--------------------------------------------------------------------------
169  Lagrangian basis function evaluation
170  --------------------------------------------------------------------------*/
171 
172 /* given the Lagrangian nodes (z,n) and evaluation points (x,m)
173  evaluate all Lagrangian basis functions at all points x
174 
175  inner index of output J is the basis function index (row-major format)
176  provide work array with space for 4*n doubles
177  */
178 void lagrange_weights( const realType* z, unsigned n, const realType* x, unsigned m, realType* J, realType* work );
179 
180 /* given the Lagrangian nodes (z,n) and evaluation points (x,m)
181  evaluate all Lagrangian basis functions and their derivatives
182 
183  inner index of outputs J,D is the basis function index (row-major format)
184  provide work array with space for 6*n doubles
185  */
186 void lagrange_weights_deriv( const realType* z,
187  unsigned n,
188  const realType* x,
189  unsigned m,
190  realType* J,
191  realType* D,
192  realType* work );
193 
194 /*--------------------------------------------------------------------------
195  Speedy Lagrangian Interpolation
196 
197  Usage:
198 
199  lagrange_data p;
200  lagrange_setup(&p,z,n); * setup for nodes z[0 ... n-1] *
201 
202  the weights
203  p->J [0 ... n-1] interpolation weights
204  p->D [0 ... n-1] 1st derivative weights
205  p->D2[0 ... n-1] 2nd derivative weights
206  are computed for a given x with:
207  lagrange_0(p,x); * compute p->J *
208  lagrange_1(p,x); * compute p->J, p->D *
209  lagrange_2(p,x); * compute p->J, p->D, p->D2 *
210  lagrange_2u(p); * compute p->D2 after call of lagrange_1(p,x); *
211  These functions use the z array supplied to setup
212  (that pointer should not be freed between calls)
213  Weights for x=z[0] and x=z[n-1] are computed during setup; access as:
214  p->J_z0, etc. and p->J_zn, etc.
215 
216  lagrange_free(&p); * deallocate memory allocated by setup *
217  --------------------------------------------------------------------------*/
218 
219 typedef struct
220 {
221  unsigned n; /* number of Lagrange nodes */
222  const realType* z; /* Lagrange nodes (user-supplied) */
223  realType *J, *D, *D2; /* weights for 0th,1st,2nd derivatives */
224  realType *J_z0, *D_z0, *D2_z0; /* ditto at z[0] (computed at setup) */
225  realType *J_zn, *D_zn, *D2_zn; /* ditto at z[n-1] (computed at setup) */
226  realType *w, *d, *u0, *v0, *u1, *v1, *u2, *v2; /* work data */
227 } lagrange_data;
228 
229 void lagrange_setup( lagrange_data* p, const realType* z, unsigned n );
230 void lagrange_free( lagrange_data* p );
231 
232 void lagrange_0( lagrange_data* p, realType x );
233 void lagrange_1( lagrange_data* p, realType x );
234 void lagrange_2( lagrange_data* p, realType x );
235 void lagrange_2u( lagrange_data* p );
236 
237 /*======================================================
238 / from tensor.h
239 /======================================================
240 */
241 
242 /*--------------------------------------------------------------------------
243  1-,2-,3-d Tensor Application
244 
245  the 3d case:
246  tensor_f3(R,mr,nr, S,ms,ns, T,mt,nt, u,v, work1,work2)
247  gives v = [ R (x) S (x) T ] u
248  where R is mr x nr, S is ms x ns, T is mt x nt,
249  each in row- or column-major format according to f := r | c
250  u is nr x ns x nt in column-major format (inner index is r)
251  v is mr x ms x mt in column-major format (inner index is r)
252  --------------------------------------------------------------------------*/
253 
254 void tensor_c1( const realType* R, unsigned mr, unsigned nr, const realType* u, realType* v );
255 void tensor_r1( const realType* R, unsigned mr, unsigned nr, const realType* u, realType* v );
256 
257 /* work holds mr*ns realTypes */
258 void tensor_c2( const realType* R,
259  unsigned mr,
260  unsigned nr,
261  const realType* S,
262  unsigned ms,
263  unsigned ns,
264  const realType* u,
265  realType* v,
266  realType* work );
267 void tensor_r2( const realType* R,
268  unsigned mr,
269  unsigned nr,
270  const realType* S,
271  unsigned ms,
272  unsigned ns,
273  const realType* u,
274  realType* v,
275  realType* work );
276 
277 /* work1 holds mr*ns*nt realTypes,
278  work2 holds mr*ms*nt realTypes */
279 void tensor_c3( const realType* R,
280  unsigned mr,
281  unsigned nr,
282  const realType* S,
283  unsigned ms,
284  unsigned ns,
285  const realType* T,
286  unsigned mt,
287  unsigned nt,
288  const realType* u,
289  realType* v,
290  realType* work1,
291  realType* work2 );
292 void tensor_r3( const realType* R,
293  unsigned mr,
294  unsigned nr,
295  const realType* S,
296  unsigned ms,
297  unsigned ns,
298  const realType* T,
299  unsigned mt,
300  unsigned nt,
301  const realType* u,
302  realType* v,
303  realType* work1,
304  realType* work2 );
305 
306 /*--------------------------------------------------------------------------
307  1-,2-,3-d Tensor Application of Row Vectors (for Interpolation)
308 
309  the 3d case:
310  v = tensor_i3(Jr,nr, Js,ns, Jt,nt, u, work)
311  same effect as tensor_r3(Jr,1,nr, Js,1,ns, Jt,1,nt, u,&v, work1,work2):
312  gives v = [ Jr (x) Js (x) Jt ] u
313  where Jr, Js, Jt are row vectors (interpolation weights)
314  u is nr x ns x nt in column-major format (inner index is r)
315  v is a scalar
316  --------------------------------------------------------------------------*/
317 
318 realType tensor_i1( const realType* Jr, unsigned nr, const realType* u );
319 
320 /* work holds ns realTypes */
321 realType tensor_i2( const realType* Jr,
322  unsigned nr,
323  const realType* Js,
324  unsigned ns,
325  const realType* u,
326  realType* work );
327 
328 /* work holds ns*nt + nt realTypes */
329 realType tensor_i3( const realType* Jr,
330  unsigned nr,
331  const realType* Js,
332  unsigned ns,
333  const realType* Jt,
334  unsigned nt,
335  const realType* u,
336  realType* work );
337 
338 /*--------------------------------------------------------------------------
339  1-,2-,3-d Tensor Application of Row Vectors
340  for simultaneous Interpolation and Gradient computation
341 
342  the 3d case:
343  v = tensor_ig3(Jr,Dr,nr, Js,Ds,ns, Jt,Dt,nt, u,g, work)
344  gives v = [ Jr (x) Js (x) Jt ] u
345  g_0 = [ Dr (x) Js (x) Jt ] u
346  g_1 = [ Jr (x) Ds (x) Jt ] u
347  g_2 = [ Jr (x) Js (x) Dt ] u
348  where Jr,Dr,Js,Ds,Jt,Dt are row vectors
349  (interpolation & derivative weights)
350  u is nr x ns x nt in column-major format (inner index is r)
351  v is a scalar, g is an array of 3 realTypes
352  --------------------------------------------------------------------------*/
353 
354 realType tensor_ig1( const realType* Jr, const realType* Dr, unsigned nr, const realType* u, realType* g );
355 
356 /* work holds 2*ns realTypes */
357 realType tensor_ig2( const realType* Jr,
358  const realType* Dr,
359  unsigned nr,
360  const realType* Js,
361  const realType* Ds,
362  unsigned ns,
363  const realType* u,
364  realType* g,
365  realType* work );
366 
367 /* work holds 2*ns*nt + 3*ns realTypes */
368 realType tensor_ig3( const realType* Jr,
369  const realType* Dr,
370  unsigned nr,
371  const realType* Js,
372  const realType* Ds,
373  unsigned ns,
374  const realType* Jt,
375  const realType* Dt,
376  unsigned nt,
377  const realType* u,
378  realType* g,
379  realType* work );
380 
381 /*======================================================
382 / from findpt.h
383 /======================================================
384 */
385 
386 typedef struct
387 {
388  realType x[2], A[4], axis_bnd[4];
389 } obbox_2;
390 
391 typedef struct
392 {
393  realType x[3], A[9], axis_bnd[6];
394 } obbox_3;
395 
396 typedef struct
397 {
398  unsigned hash_n;
399  realType bnd[4]; /* bounds for all elements */
400  realType fac[2]; /* fac[i] = hash_n / (bnd[2*i+1]-bnd[2*i]) */
401  obbox_2* obb; /* obb[nel] -- bounding box info for each element */
402  uint* offset; /* hash table -- for cell i,j:
403  uint index = j*hash_n+i,
404  b = offset[index ],
405  e = offset[index+1];
406  elements in cell are
407  offset[b], offset[b+1], ..., offset[e-1] */
408  unsigned max; /* maximum # of elements in any cell */
409 } hash_data_2;
410 
411 typedef struct
412 {
413  unsigned hash_n;
414  realType bnd[6]; /* bounds for all elements */
415  realType fac[3]; /* fac[i] = hash_n / (bnd[2*i+1]-bnd[2*i]) */
416  obbox_3* obb; /* obb[nel] -- bounding box info for each element */
417  uint* offset; /* hash table -- for cell i,j,k:
418  uint index = (k*hash_n+j)*hash_n+i,
419  b = offset[index ],
420  e = offset[index+1];
421  elements in cell are
422  offset[b], offset[b+1], ..., offset[e-1] */
423  unsigned max; /* maximum # of elements in any cell */
424 } hash_data_3;
425 
426 typedef struct
427 {
428  uint el;
429  realType r[3];
430  realType dist;
431 } findpt_listel;
432 
433 typedef struct
434 {
435  unsigned constraints;
436  unsigned de, d1;
437  realType *x[2], *fd1[2];
438 } opt_edge_data_2;
439 
440 typedef struct
441 {
442  unsigned constraints;
443  realType x[2], jac[4];
444 } opt_point_data_2;
445 
446 typedef struct
447 {
448  lagrange_data* ld;
449  unsigned size[3];
450  const realType* elx[2];
451  opt_edge_data_2 ed;
452  opt_point_data_2 pd;
453  realType* work;
454  realType x[2], jac[4];
455 } opt_data_2;
456 
457 typedef struct
458 {
459  const realType* xw[2]; /* geometry data */
460  realType* z[2]; /* lobatto nodes */
461  lagrange_data ld[2]; /* interpolation, derivative weights & data */
462  unsigned nptel; /* nodes per element */
463  hash_data_2* hash; /* geometric hashing data */
464  findpt_listel *list, **sorted, **end;
465  /* pre-allocated list of elements to
466  check (found by hashing), and
467  pre-allocated list of pointers into
468  the first list for sorting */
469  opt_data_2* od; /* data for the optimization algorithm */
470  realType* od_work;
471 } findpt_data_2;
472 
473 typedef struct
474 {
475  unsigned constraints;
476  unsigned dn, d1, d2;
477  realType *x[3], *fdn[3];
478 } opt_face_data_3;
479 
480 typedef struct
481 {
482  unsigned constraints;
483  unsigned de, d1, d2;
484  realType *x[3], *fd1[3], *fd2[3];
485 } opt_edge_data_3;
486 
487 typedef struct
488 {
489  unsigned constraints;
490  realType x[3], jac[9];
491 } opt_point_data_3;
492 
493 typedef struct
494 {
495  lagrange_data* ld;
496  unsigned size[4];
497  const realType* elx[3];
498  opt_face_data_3 fd;
499  opt_edge_data_3 ed;
500  opt_point_data_3 pd;
501  realType* work;
502  realType x[3], jac[9];
503 } opt_data_3;
504 
505 typedef struct
506 {
507  const realType* xw[3]; /* geometry data */
508  realType* z[3]; /* lobatto nodes */
509  lagrange_data ld[3]; /* interpolation, derivative weights & data */
510  unsigned nptel; /* nodes per element */
511  hash_data_3* hash; /* geometric hashing data */
512  findpt_listel *list, **sorted, **end;
513  /* pre-allocated list of elements to
514  check (found by hashing), and
515  pre-allocated list of pointers into
516  the first list for sorting */
517  opt_data_3* od; /* data for the optimization algorithm */
518  realType* od_work;
519 } findpt_data_3;
520 
521 findpt_data_2* findpt_setup_2( const realType* const xw[2],
522  const unsigned n[2],
523  uint nel,
524  uint max_hash_size,
525  realType bbox_tol );
526 findpt_data_3* findpt_setup_3( const realType* const xw[3],
527  const unsigned n[3],
528  uint nel,
529  uint max_hash_size,
530  realType bbox_tol );
531 
532 void findpt_free_2( findpt_data_2* p );
533 void findpt_free_3( findpt_data_3* p );
534 
535 const realType* findpt_allbnd_2( const findpt_data_2* p );
536 const realType* findpt_allbnd_3( const findpt_data_3* p );
537 
538 typedef int ( *findpt_func )( void*, const realType*, int, uint*, realType*, realType* );
539 int findpt_2( findpt_data_2* p, const realType x[2], int guess, uint* el, realType r[2], realType* dist );
540 int findpt_3( findpt_data_3* p, const realType x[3], int guess, uint* el, realType r[3], realType* dist );
541 
542 void findpt_weights_2( findpt_data_2* p, const realType r[2] );
543 
544 void findpt_weights_3( findpt_data_3* p, const realType r[3] );
545 
546 double findpt_eval_2( findpt_data_2* p, const realType* u );
547 
548 double findpt_eval_3( findpt_data_3* p, const realType* u );
549 
550 /*======================================================
551 / from extrafindpt.h
552 /======================================================
553 */
554 
555 void opt_alloc_3( opt_data_3* p, lagrange_data* ld );
556 void opt_free_3( opt_data_3* p );
557 double opt_findpt_3( opt_data_3* p,
558  const realType* const elx[3],
559  const realType xstar[3],
560  realType r[3],
561  unsigned* constr );
562 void opt_vol_set_intp_3( opt_data_3* p, const realType r[3] );
563 
564 /* for 2d spectralQuad */
565 /*--------------------------------------------------------------------------
566 
567  2 - D
568 
569  --------------------------------------------------------------------------*/
570 
571 void opt_alloc_2( opt_data_2* p, lagrange_data* ld );
572 void opt_free_2( opt_data_2* p );
573 double opt_findpt_2( opt_data_2* p,
574  const realType* const elx[2],
575  const realType xstar[2],
576  realType r[2],
577  unsigned* constr );
578 
579 extern const unsigned opt_no_constraints_2;
580 extern const unsigned opt_no_constraints_3;
581 
582 /*======================================================
583 / from errmem.h
584 /======================================================
585 */
586 
587 /* requires:
588  <stdlib.h> for malloc, calloc, realloc, free
589 */
590 
591 /*--------------------------------------------------------------------------
592  Error Reporting
593  Memory Allocation Wrappers to Catch Out-of-memory
594  --------------------------------------------------------------------------*/
595 
596 /* #include "malloc.h" */
597 #include <stdlib.h>
598 
599 #ifdef __GNUC__
600 void fail( const char* fmt, ... ) __attribute__( ( noreturn ) );
601 #define MAYBE_UNUSED __attribute__( ( unused ) )
602 #else
603 void fail( const char* fmt, ... );
604 #define MAYBE_UNUSED
605 #endif
606 
607 #if 0
608 {}
609 #endif
610 
611 static void* smalloc( size_t size, const char* file ) MAYBE_UNUSED;
612 static void* smalloc( size_t size, const char* file )
613 {
614  void* res = malloc( size );
615  if( !res && size ) fail( "%s: allocation of %d bytes failed\n", file, (int)size );
616  return res;
617 }
618 
619 static void* scalloc( size_t nmemb, size_t size, const char* file ) MAYBE_UNUSED;
620 static void* scalloc( size_t nmemb, size_t size, const char* file )
621 {
622  void* res = calloc( nmemb, size );
623  if( !res && nmemb ) fail( "%s: allocation of %d bytes failed\n", file, (int)size * nmemb );
624  return res;
625 }
626 
627 static void* srealloc( void* ptr, size_t size, const char* file ) MAYBE_UNUSED;
628 static void* srealloc( void* ptr, size_t size, const char* file )
629 {
630  void* res = realloc( ptr, size );
631  if( !res && size ) fail( "%s: allocation of %d bytes failed\n", file, (int)size );
632  return res;
633 }
634 
635 #define tmalloc( type, count ) ( (type*)smalloc( ( count ) * sizeof( type ), __FILE__ ) )
636 #define tcalloc( type, count ) ( (type*)scalloc( ( count ), sizeof( type ), __FILE__ ) )
637 #define trealloc( type, ptr, count ) ( (type*)srealloc( ( ptr ), ( count ) * sizeof( type ), __FILE__ ) )
638 /*
639 typedef struct { size_t size; void *ptr; } buffer;
640 
641 static void buffer_init_(buffer *b, size_t size, const char *file) MAYBE_UNUSED;
642 static void buffer_init_(buffer *b, size_t size, const char *file)
643 {
644  b->size=size, b->ptr=smalloc(size,file);
645 }
646 static void buffer_reserve_(buffer *b, size_t min, const char *file)
647  MAYBE_UNUSED;
648 static void buffer_reserve_(buffer *b, size_t min, const char *file)
649 {
650  size_t size = b->size;
651  if(size<min) {
652  size+=size/2+1;
653  if(size<min) size=min;
654  b->ptr=srealloc(b->ptr,size,file);
655  }
656 }
657 static void buffer_free(buffer *b) MAYBE_UNUSED;
658 static void buffer_free(buffer *b) { free(b->ptr); }
659 
660 #define buffer_init(b,size) buffer_init_(b,size,__FILE__)
661 #define buffer_reserve(b,min) buffer_reserve_(b,min,__FILE__)
662 */
663 
664 /*======================================================
665 / from minmax.h
666 /======================================================
667 */
668 
669 /*--------------------------------------------------------------------------
670  Min, Max, Norm
671  --------------------------------------------------------------------------*/
672 
673 #ifdef __GNUC__
674 #define MAYBE_UNUSED __attribute__( ( unused ) )
675 #else
676 #define MAYBE_UNUSED
677 #endif
678 
679 #define DECLMINMAX( type, prefix ) \
680  static type prefix##min_2( type a, type b ) MAYBE_UNUSED; \
681  static type prefix##min_2( type a, type b ) \
682  { \
683  return b < a ? b : a; \
684  } \
685  static type prefix##max_2( type a, type b ) MAYBE_UNUSED; \
686  static type prefix##max_2( type a, type b ) \
687  { \
688  return a > b ? a : b; \
689  } \
690  static void prefix##minmax_2( type* min, type* max, type a, type b ) MAYBE_UNUSED; \
691  static void prefix##minmax_2( type* min, type* max, type a, type b ) \
692  { \
693  if( b < a ) \
694  *min = b, *max = a; \
695  else \
696  *min = a, *max = b; \
697  } \
698  static type prefix##min_3( type a, type b, type c ) MAYBE_UNUSED; \
699  static type prefix##min_3( type a, type b, type c ) \
700  { \
701  return b < a ? ( c < b ? c : b ) : ( c < a ? c : a ); \
702  } \
703  static type prefix##max_3( type a, type b, type c ) MAYBE_UNUSED; \
704  static type prefix##max_3( type a, type b, type c ) \
705  { \
706  return a > b ? ( a > c ? a : c ) : ( b > c ? b : c ); \
707  } \
708  static void prefix##minmax_3( type* min, type* max, type a, type b, type c ) MAYBE_UNUSED; \
709  static void prefix##minmax_3( type* min, type* max, type a, type b, type c ) \
710  { \
711  if( b < a ) \
712  *min = prefix##min_2( b, c ), *max = prefix##max_2( a, c ); \
713  else \
714  *min = prefix##min_2( a, c ), *max = prefix##max_2( b, c ); \
715  }
716 
717 DECLMINMAX( int, i )
718 DECLMINMAX( unsigned, u )
719 DECLMINMAX( realType, r )
720 #undef DECLMINMAX
721 
722 static realType r1norm_1( realType a ) MAYBE_UNUSED;
723 static realType r1norm_1( realType a )
724 {
725  return mbabs( a );
726 }
727 static realType r1norm_2( realType a, realType b ) MAYBE_UNUSED;
728 static realType r1norm_2( realType a, realType b )
729 {
730  return mbabs( a ) + mbabs( b );
731 }
732 static realType r1norm_3( realType a, realType b, realType c ) MAYBE_UNUSED;
733 static realType r1norm_3( realType a, realType b, realType c )
734 {
735  return mbabs( a ) + mbabs( b ) + mbabs( c );
736 }
737 static realType r2norm_1( realType a ) MAYBE_UNUSED;
738 static realType r2norm_1( realType a )
739 {
740  return mbsqrt( a * a );
741 }
742 static realType r2norm_2( realType a, realType b ) MAYBE_UNUSED;
743 static realType r2norm_2( realType a, realType b )
744 {
745  return mbsqrt( a * a + b * b );
746 }
747 static realType r2norm_3( realType a, realType b, realType c ) MAYBE_UNUSED;
748 static realType r2norm_3( realType a, realType b, realType c )
749 {
750  return mbsqrt( a * a + b * b + c * c );
751 }
752 
753 #endif