HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
arrayfunctions.h
Go to the documentation of this file.
1 
6 #ifndef _ARRAYFUNCTIONS_H_
7 #define _ARRAYFUNCTIONS_H_
8 
9 #include <stdio.h>
10 #include <stdlib.h>
11 #include <basic.h>
12 #include <math.h>
13 
21 #define _ArrayIndexnD_(N,index,imax,i,ghost) \
22  { \
23  int arraycounter, term=1, index_copy=index; \
24  for (arraycounter=0; arraycounter<N; arraycounter++) term *= (imax[arraycounter]+2*(ghost));\
25  for (arraycounter=N-1; arraycounter>=0; arraycounter--) {\
26  term /= (imax[arraycounter]+2*(ghost)); \
27  i[arraycounter] = index_copy/term; \
28  index_copy -= i[arraycounter]*term; \
29  } \
30  for (arraycounter=0; arraycounter<N; arraycounter++) i[arraycounter] -= (ghost);\
31  }
32 
40 #define _ArrayIndex1D_(N,imax,i,ghost,index) \
41  { \
42  index = i[N-1]+(ghost); \
43  int arraycounter; \
44  for (arraycounter = (N)-2; arraycounter > -1; arraycounter--) { \
45  index = ((index*(imax[arraycounter]+2*(ghost))) + (i[arraycounter]+(ghost))); \
46  } \
47  }
48 
56 #define _ArrayIndex1D2_(N,imax,i,ghost,index) \
57  { \
58  index = i[1]+(ghost); \
59  index = ((index*(imax[0]+2*(ghost))) + (i[0]+(ghost))); \
60  }
61 
69 #define _ArrayIndex1D3_(N,imax,i,ghost,index) \
70  { \
71  index = i[2]+(ghost); \
72  index = ((index*(imax[1]+2*(ghost))) + (i[1]+(ghost))); \
73  index = ((index*(imax[0]+2*(ghost))) + (i[0]+(ghost))); \
74  }
75 
83 #define _ArrayIndex1DWO_(N,imax,i,offset,ghost,index) \
84  { \
85  index = i[N-1]+(ghost)+ offset[N-1];\
86  int arraycounter; \
87  for (arraycounter = (N)-2; arraycounter > -1; arraycounter--) { \
88  index = ((index*(imax[arraycounter]+2*(ghost))) + (i[arraycounter]+(ghost)+offset[arraycounter]));\
89  } \
90  }
91 
100 #define _ArrayIndex1DWO2_(N,imax,i,offset,ghost,index) \
101  { \
102  index = i[1]+(ghost)+ offset[1];\
103  index = ((index*(imax[0]+2*(ghost))) + (i[0]+(ghost)+offset[0]));\
104  }
105 
114 #define _ArrayIndex1DWO3_(N,imax,i,offset,ghost,index) \
115  { \
116  index = i[2]+(ghost)+ offset[2];\
117  index = ((index*(imax[1]+2*(ghost))) + (i[1]+(ghost)+offset[1]));\
118  index = ((index*(imax[0]+2*(ghost))) + (i[0]+(ghost)+offset[0]));\
119  }
120 
126 #define _ArrayIncrementIndex_(N,imax,i,done) \
127  { \
128  int arraycounter = 0; \
129  while (arraycounter < (N)) { \
130  if (i[arraycounter] == imax[arraycounter]-1) { \
131  i[arraycounter] = 0; \
132  arraycounter++; \
133  } else { \
134  i[arraycounter]++; \
135  break; \
136  } \
137  } \
138  if (arraycounter == (N)) done = 1; \
139  else done = 0; \
140  }
141 
149 #define _ArrayIncrementIndexWithLBound_(N,imax,imin,i,done) \
150  { \
151  int arraycounter = 0; \
152  while (arraycounter < (N)) { \
153  if (i[arraycounter] == imax[arraycounter]-1) { \
154  i[arraycounter] = imin[arraycounter]; \
155  arraycounter++; \
156  } else { \
157  i[arraycounter]++; \
158  break; \
159  } \
160  } \
161  if (arraycounter == (N)) done = 1; \
162  else done = 0; \
163  }
164 
169 #define _ArraySetValue_(x,size,value) \
170  { \
171  int arraycounter; \
172  for (arraycounter = 0; arraycounter < (size); arraycounter++) x[arraycounter] = (value); \
173  }
174 
179 #define _ArraySum1D_(x,a,size) \
180  { \
181  a = 0; \
182  int arraycounter; \
183  for (arraycounter=0; arraycounter<size; arraycounter++) a += x[arraycounter]; \
184  }
185 
190 #define _ArrayScale1D_(x,a,size) \
191  { \
192  int arraycounter; \
193  for (arraycounter=0; arraycounter<size; arraycounter++) x[arraycounter] *= a; \
194  }
195 
201 #define _ArraySubtract1D_(x,a,b,size) \
202  { \
203  int arraycounter; \
204  for (arraycounter=0; arraycounter<size; arraycounter++) x[arraycounter] = a[arraycounter] - b[arraycounter]; \
205  }
206 
212 #define _ArrayAdd1D_(x,a,b,size) \
213  { \
214  int arraycounter; \
215  for (arraycounter=0; arraycounter<size; arraycounter++) x[arraycounter] = a[arraycounter] + b[arraycounter]; \
216  }
217 
223 #define _ArrayMultiply1D_(x,a,b,size) \
224  { \
225  int arraycounter; \
226  for (arraycounter=0; arraycounter<size; arraycounter++) x[arraycounter] = a[arraycounter] * b[arraycounter]; \
227  }
228 
234 #define _ArrayMultiply3Add1D_(x,a,b,c,d,e,f,size) \
235  { \
236  int arraycounter; \
237  for (arraycounter=0; arraycounter<size; arraycounter++) \
238  x[arraycounter] = a[arraycounter]*b[arraycounter]+c[arraycounter]*d[arraycounter]+e[arraycounter]*f[arraycounter]; \
239  }
240 
245 #define _ArrayConvexCombination1D_(z,a,x,y,size) \
246  { \
247  int arraycounter; \
248  for (arraycounter=0; arraycounter<size; arraycounter++) \
249  z[arraycounter] = a[arraycounter]*x[arraycounter]+(1.0-a[arraycounter])*y[arraycounter]; \
250  }
251 
258 #define _ArrayAYPX_(x,a,y,size) \
259  { \
260  int arraycounter; \
261  for (arraycounter=0; arraycounter<size; arraycounter++) \
262  y[arraycounter] = a*y[arraycounter] + x[arraycounter];\
263  }
264 
271 #define _ArrayAXPY_(x,a,y,size) \
272  { \
273  int arraycounter; \
274  for (arraycounter=0; arraycounter<size; arraycounter++) y[arraycounter] += a*x[arraycounter]; \
275  }
276 
283 #define _ArrayAXBY_(z,a,x,b,y,size) \
284  { \
285  int arraycounter; \
286  for (arraycounter=0; arraycounter<size; arraycounter++) z[arraycounter] = a*x[arraycounter]+b*y[arraycounter]; \
287  }
288 
295 #define _ArrayAXBYCZ_(w,a,x,b,y,c,z,size) \
296  { \
297  int arraycounter; \
298  for (arraycounter=0; arraycounter<size; arraycounter++) w[arraycounter] = a*x[arraycounter]+b*y[arraycounter]+c*z[arraycounter]; \
299  }
300 
307 #define _ArrayScaledAXPY_(x,a,e,y,size) \
308  { \
309  int arraycounter; \
310  for (arraycounter=0; arraycounter<size; arraycounter++) \
311  y[arraycounter] = e*(y[arraycounter]+a*x[arraycounter]); \
312  }
313 
319 #define _ArrayCopy1D_(x,y,size) \
320  { \
321  int arraycounter; \
322  for (arraycounter = 0; arraycounter < size; arraycounter++) y[arraycounter] = x[arraycounter]; \
323  }
324 
331 #define _ArrayCopy1D2_(x,y,size) \
332  { \
333  y[0] = x[0]; \
334  y[1] = x[1]; \
335  }
336 
343 #define _ArrayCopy1D3_(x,y,size) \
344  { \
345  y[0] = x[0]; \
346  y[1] = x[1]; \
347  y[2] = x[2]; \
348  }
349 
356 #define _ArrayScaleCopy1D_(x,a,y,size) \
357  { \
358  int arraycounter; \
359  for (arraycounter = 0; arraycounter < size; arraycounter++) y[arraycounter] = a*x[arraycounter]; \
360  }
361 
368 #define _ArrayAddCopy1D_(x,a,y,size) \
369  { \
370  int arraycounter; \
371  for (arraycounter = 0; arraycounter < size; arraycounter++) y[arraycounter] = a+x[arraycounter]; \
372  }
373 
377 #define _ArrayProduct1D_(x,size,p) \
378  { \
379  int arraycounter = 0; p = 1; \
380  for (arraycounter=0; arraycounter<size; arraycounter++) p *= x[arraycounter]; \
381  }
382 
388 #define _ArrayBlockMultiply_(x,a,n,bs) \
389  { \
390  int arraycounter1,arraycounter2; \
391  for (arraycounter1=0; arraycounter1<n; arraycounter1++) { \
392  for (arraycounter2=0; arraycounter2<bs; arraycounter2++) { \
393  x[bs*arraycounter1+arraycounter2] *= a[arraycounter1]; \
394  } \
395  }\
396  }
397 
398 #if !defined(INLINE)
399 # define INLINE inline
400 #endif
401 
402 INLINE int ArrayCopynD (int,const double*,double*,int*,int,int,int*,int);
403 INLINE double ArrayMaxnD (int,int,int*,int,int*,double*);
404 INLINE double ArraySumSquarenD (int,int,int*,int,int*,double*);
405 INLINE double ArraySumAbsnD (int,int,int*,int,int*,double*);
406 INLINE void ArrayCheckEqual (const char *,const double *,const double *,int,int);
407 INLINE void ArrayCheckEqual2 (const char *,const double *,const double *,int);
408 
410 INLINE int ArrayCopynD(int ndims,
411  const double *x,
412  double *y,
413  int *dim,
414  int g1,
415  int g2,
416  int *index,
417  int nvars
419  )
420 {
421  if (!y) {
422  fprintf(stderr,"Error in ArrayCopynD(): array \"y\" not allocated.\n");
423  return(1);
424  }
425  if (!x) {
426  fprintf(stderr,"Error in ArrayCopynD(): array \"x\" not allocated.\n");
427  return(1);
428  }
429  int done = 0;
430  _ArraySetValue_(index,ndims,0);
431  while (!done) {
432  int p1, p2;
433  _ArrayIndex1D_(ndims,dim,index,g1,p1);
434  _ArrayIndex1D_(ndims,dim,index,g2,p2);
435  _ArrayCopy1D_((x+p1*nvars),(y+p2*nvars),nvars);
436  _ArrayIncrementIndex_(ndims,dim,index,done);
437  }
438  return(0);
439 }
440 
444 INLINE int ArrayCopynDComponent( int ndims,
445  const double *x,
446  double *y,
447  int *dim,
448  int g1,
449  int g2,
450  int *index,
451  int nvars_from,
452  int nvars_to,
453  int var_from,
454  int var_to )
455 {
456  if (!y) {
457  fprintf(stderr,"Error in ArrayCopynD(): array \"y\" not allocated.\n");
458  return(1);
459  }
460  if (!x) {
461  fprintf(stderr,"Error in ArrayCopynD(): array \"x\" not allocated.\n");
462  return(1);
463  }
464  if ((var_from >= nvars_from) || (var_from < 0)) {
465  fprintf(stderr,"Error in ArrayCopynD(): specified component exceeds bounds.\n");
466  return(1);
467  }
468  if ((var_to >= nvars_to) || (var_to < 0)) {
469  fprintf(stderr,"Error in ArrayCopynD(): specified component exceeds bounds.\n");
470  return(1);
471  }
472 
473  int done = 0;
474  _ArraySetValue_(index,ndims,0);
475  while (!done) {
476  int p1, p2;
477  _ArrayIndex1D_(ndims,dim,index,g1,p1);
478  _ArrayIndex1D_(ndims,dim,index,g2,p2);
479  (y+p2*nvars_to)[var_to] = (x+p1*nvars_from)[var_from];
480  _ArrayIncrementIndex_(ndims,dim,index,done);
481  }
482  return(0);
483 }
484 
486 INLINE double ArrayMaxnD(int nvars,
488  int ndims,
489  int *dim,
490  int ghosts,
491  int *index,
492  double *x
493  )
494 {
495  double sum = 0;
496  int done = 0; _ArraySetValue_(index,ndims,0);
497  while (!done) {
498  int p; _ArrayIndex1D_(ndims,dim,index,ghosts,p);
499  int v;
500  for (v=0; v<nvars; v++) {
501  double term = ( x[p*nvars+v]>0 ? x[p*nvars+v] : -x[p*nvars+v] );
502  if (term > sum) sum = term;
503  }
504  _ArrayIncrementIndex_(ndims,dim,index,done);
505  }
506  return(sum);
507 }
508 
510 INLINE double ArraySumAbsnD(int nvars,
511  int ndims,
512  int *dim,
513  int ghosts,
514  int *index,
515  double *x
516  )
517 {
518  double sum = 0;
519  int done = 0; _ArraySetValue_(index,ndims,0);
520  while (!done) {
521  int p; _ArrayIndex1D_(ndims,dim,index,ghosts,p);
522  int v; for (v=0; v<nvars; v++) sum += ( x[p*nvars+v]>0 ? x[p*nvars+v] : -x[p*nvars+v] );
523  _ArrayIncrementIndex_(ndims,dim,index,done);
524  }
525  return(sum);
526 }
527 
529 INLINE double ArraySumSquarenD(int nvars,
530  int ndims,
531  int *dim,
532  int ghosts,
533  int *index,
534  double *x
535  )
536 {
537  double sum = 0;
538  int done = 0; _ArraySetValue_(index,ndims,0);
539  while (!done) {
540  int p; _ArrayIndex1D_(ndims,dim,index,ghosts,p);
541  int v; for (v=0; v<nvars; v++) sum += (x[p*nvars+v]*x[p*nvars+v]);
542  _ArrayIncrementIndex_(ndims,dim,index,done);
543  }
544  return(sum);
545 }
546 
547 INLINE void ArrayCheckEqual(const char *msg, const double *var_adj, const double *var_sep, int npoints, int nvars)
548 {
549  double max_err = 0.0;
550  for (int j=0; j<npoints; j++) {
551  for (int v=0; v<nvars; v++) {
552  if (var_sep[j+v*npoints] != var_adj[j*nvars+v]) {
553  double diff = fabs(var_sep[j+v*npoints]-var_adj[j*nvars+v]);
554  max_err = (max_err < diff) ? diff : max_err;
555  }
556  }
557  }
558 
559  if (max_err > 1e-10) {
560  printf("ArrayCheckEqual: %-30s max_err = %e\n", msg, max_err);
561  exit(0);
562  }
563  return;
564 }
565 
566 INLINE void ArrayCheckEqual2(const char *msg, const double *var_adj, const double *var_sep, int npoints)
567 {
568  double max_err = 0.0;
569  for (int j=0; j<npoints; j++) {
570  if (var_sep[j] != var_adj[j]) {
571  double diff = fabs(var_sep[j]-var_adj[j]);
572  max_err = (max_err < diff) ? diff : max_err;
573  }
574  }
575 
576  if (max_err > 1e-10) {
577  printf("ArrayCheckEqual: %-30s max_err = %e\n", msg, max_err);
578  exit(0);
579  }
580  return;
581 }
582 
583 
584 #endif
StdVecOps::sum
long sum(const std::vector< int > &a_iv)
Definition: std_vec_ops.h:18
basic.h
Some basic definitions and macros.
ArrayCheckEqual
INLINE void ArrayCheckEqual(const char *, const double *, const double *, int, int)
Definition: arrayfunctions.h:547
ArraySumAbsnD
INLINE double ArraySumAbsnD(int, int, int *, int, int *, double *)
Definition: arrayfunctions.h:510
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
ArrayCopynDComponent
INLINE int ArrayCopynDComponent(int ndims, const double *x, double *y, int *dim, int g1, int g2, int *index, int nvars_from, int nvars_to, int var_from, int var_to)
Definition: arrayfunctions.h:444
ArrayCopynD
INLINE int ArrayCopynD(int, const double *, double *, int *, int, int, int *, int)
Definition: arrayfunctions.h:410
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
ArrayMaxnD
INLINE double ArrayMaxnD(int, int, int *, int, int *, double *)
Definition: arrayfunctions.h:486
ArrayCheckEqual2
INLINE void ArrayCheckEqual2(const char *, const double *, const double *, int)
Definition: arrayfunctions.h:566
INLINE
#define INLINE
Definition: arrayfunctions.h:399
ArraySumSquarenD
INLINE double ArraySumSquarenD(int, int, int *, int, int *, double *)
Definition: arrayfunctions.h:529
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319