 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
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
1.9.1.