Nektar++
Vmath.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: Vmath.hpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Collection of templated functions for vector mathematics
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <LibUtilities/BasicUtils/Vmath.hpp>
36 #include <LibUtilities/BasicConst/NektarUnivTypeDefs.hpp>
37 #include <LibUtilities/LibUtilitiesDeclspec.h>
38 
39 namespace Vmath
40 {
41 
42  /***************** Math routines ***************/
43 
44  /// \brief Fill a vector with a constant value
45  template<class T> void Fill( int n, const T alpha, T *x, const int incx )
46  {
47  while( n-- )
48  {
49  *x = alpha;
50  x += incx;
51  }
52  }
53 
54  template LIB_UTILITIES_EXPORT void Fill( int n, const Nektar::NekDouble alpha, Nektar::NekDouble *x, const int incx );
55 
56  #define IM1 2147483563
57  #define IM2 2147483399
58  #define AM (1.0/IM1)
59  #define IMM1 (IM1-1)
60  #define IA1 40014
61  #define IA2 40692
62  #define IQ1 53668
63  #define IQ2 52774
64  #define IR1 12211
65  #define IR2 3791
66  #define NTAB 32
67  #define NDIV (1+IMM1/NTAB)
68  #define EPS 1.2e-7
69  #define RNMX (1.0-EPS)
70 
71  /// \brief Generates a number from ~Normal(0,1)
72  template<class T> T ran2 (long* idum)
73  /* ------------------------------------------------------------------------- *
74  * Ran2 from NR 2e. Returns a uniform random deviate between 0.0 &
75  * 1.0 (exclusive of endpoints). Call with idum a negative integer to
76  * initialize; thereafter, do not alter idum between successive
77  * deviates in a sequence. RNMX should approximate the largest
78  * floating value that is less than 1.
79  * ------------------------------------------------------------------------- */
80  {
81  int j;
82  long k;
83  static long idum2=123456789;
84  static long iy=0;
85  static long iv[NTAB];
86  T temp;
87 
88  if (*idum <= 0) {
89  if (-(*idum) < 1) *idum = 1;
90  else *idum = -(*idum);
91  idum2 = (*idum);
92  for (j=NTAB+7; j>=0; j--) {
93  k = (*idum) / IQ1;
94  *idum = IA1 * (*idum - k*IQ1) - k*IR1;
95  if (*idum < 0) *idum += IM1;
96  if (j < NTAB) iv[j] = *idum;
97  }
98  iy = iv[0];
99  }
100 
101  k = (*idum) / IQ1;
102  *idum = IA1*(*idum - k*IQ1) - k*IR1;
103  if (*idum < 0) *idum += IM1;
104 
105  k = idum2 / IQ2;
106  idum2 = IA2*(idum2 - k*IQ2) - k*IR2;
107  if (idum2 < 0) idum2 += IM2;
108 
109  j = iy / NDIV;
110  iy = iv[j] - idum2;
111  iv[j] = *idum;
112  if (iy < 1) iy += IMM1;
113 
114  if ((temp=AM*iy) > RNMX) return RNMX;
115  else return temp;
116  }
117 
118  #undef IM1
119  #undef IM2
120  #undef AM
121  #undef IMM1
122  #undef IA1
123  #undef IA2
124  #undef IQ1
125  #undef IQ2
126  #undef IR1
127  #undef IR2
128  #undef NTAB
129  #undef NDIV
130  #undef EPS
131  #undef RNMX
132 
133 #ifdef NEKTAR_USE_THREAD_SAFETY
134  static boost::mutex mutex;
135 #endif
136  template LIB_UTILITIES_EXPORT Nektar::NekDouble ran2 (long* idum);
137 
138  /// \brief Fills a vector with white noise.
139  template<class T> void FillWhiteNoise( int n, const T eps, T *x,
140  const int incx, int outseed)
141  {
142 #ifdef NEKTAR_USE_THREAD_SAFETY
143  // Protect the static vars here and in ran2
144  boost::mutex::scoped_lock l(mutex);
145 #endif
146 
147  // Define static variables for generating random numbers
148  static int iset = 0;
149  static T gset;
150  static long seed = 0;
151 
152  // Bypass seed if outseed was specified
153  if( outseed != 9999)
154  {
155  seed = long(outseed);
156  }
157 
158  while( n-- )
159  {
160  T fac, rsq, v1, v2;
161 
162  if (iset == 0)
163  {
164  do
165  {
166  v1 = 2.0 * ran2<T> (&seed) - 1.0;
167  v2 = 2.0 * ran2<T> (&seed) - 1.0;
168  rsq = v1*v1 + v2*v2;
169  } while (rsq >= 1.0 || rsq == 0.0);
170  fac = sqrt(-2.0 * log (rsq) / rsq);
171  gset = v1 * fac;
172  iset = 1;
173  *x = eps * v2 * fac;
174  }
175  else
176  {
177  iset = 0;
178  *x = eps * gset;
179  }
180  x += incx;
181  }
182  }
183  template LIB_UTILITIES_EXPORT void FillWhiteNoise( int n, const Nektar::NekDouble eps, Nektar::NekDouble *x, const int incx, int outseed);
184 
185  /// \brief Multiply vector z = x*y
186  template<class T> void Vmul( int n, const T *x, const int incx, const T *y,
187  const int incy, T*z, const int incz)
188  {
189  ++n;
190  if (incx == 1 && incy == 1 && incz == 1)
191  {
192  while( --n )
193  {
194  *z = (*x) * (*y);
195  ++x;
196  ++y;
197  ++z;
198  }
199  }
200  else
201  {
202  while( --n )
203  {
204  *z = (*x) * (*y);
205  x += incx;
206  y += incy;
207  z += incz;
208  }
209  }
210  }
211  template LIB_UTILITIES_EXPORT void Vmul( int n, const Nektar::NekDouble *x, const int incx, const Nektar::NekDouble *y,
212  const int incy, Nektar::NekDouble*z, const int incz);
213 
214  /// \brief Scalar multiply y = alpha*y
215 
216  template<class T> void Smul( int n, const T alpha, const T *x, const int incx,
217  T *y, const int incy)
218  {
219  ++n;
220  if (incx == 1 && incy == 1)
221  {
222  while( --n )
223  {
224  *y = alpha * (*x);
225  ++x;
226  ++y;
227  }
228  }
229  else
230  {
231  while( --n )
232  {
233  *y = alpha * (*x);
234  x += incx;
235  y += incy;
236  }
237  }
238  }
239 
240  template LIB_UTILITIES_EXPORT void Smul( int n, const Nektar::NekDouble alpha, const Nektar::NekDouble *x, const int incx,
241  Nektar::NekDouble *y, const int incy);
242 
243  /// \brief Multiply vector z = x/y
244  template<class T> void Vdiv( int n, const T *x, const int incx, const T *y,
245  const int incy, T*z, const int incz)
246  {
247  ++n;
248  if (incx == 1 && incy == 1)
249  {
250  while( --n )
251  {
252  *z = (*x) / (*y);
253  ++x;
254  ++y;
255  ++z;
256  }
257  }
258  else
259  {
260  while( --n )
261  {
262  *z = (*x) / (*y);
263  x += incx;
264  y += incy;
265  z += incz;
266  }
267  }
268  }
269 
270  template LIB_UTILITIES_EXPORT void Vdiv( int n, const Nektar::NekDouble *x, const int incx, const Nektar::NekDouble *y,
271  const int incy, Nektar::NekDouble*z, const int incz);
272 
273  /// \brief Scalar multiply y = alpha/y
274  template<class T> void Sdiv( int n, const T alpha, const T *x,
275  const int incx, T *y, const int incy)
276  {
277  ++n;
278  if (incx == 1 && incy == 1)
279  {
280  while( --n )
281  {
282  *y = alpha / (*x);
283  ++x;
284  ++y;
285  }
286  }
287  else
288  {
289  while( --n )
290  {
291  *y = alpha / (*x);
292  x += incx;
293  y += incy;
294  }
295  }
296  }
297 
298  template LIB_UTILITIES_EXPORT void Sdiv( int n, const Nektar::NekDouble alpha, const Nektar::NekDouble *x,
299  const int incx, Nektar::NekDouble *y, const int incy);
300 
301  /// \brief Add vector z = x+y
302  template<class T> void Vadd( int n, const T *x, const int incx, const T *y,
303  const int incy, T *z, const int incz)
304  {
305  while( n-- )
306  {
307  *z = (*x) + (*y);
308  x += incx;
309  y += incy;
310  z += incz;
311  }
312  }
313 
314  template LIB_UTILITIES_EXPORT void Vadd( int n, const Nektar::NekDouble *x, const int incx, const Nektar::NekDouble *y,
315  const int incy, Nektar::NekDouble *z, const int incz);
316 
317  /// \brief Add vector y = alpha + x
318  template<class T> void Sadd( int n, const T alpha, const T *x,
319  const int incx, T *y, const int incy)
320  {
321  ++n;
322  if (incx == 1 && incy == 1)
323  {
324  while( --n )
325  {
326  *y = alpha + (*x);
327  ++x;
328  ++y;
329  }
330  }
331  else
332  {
333  while( --n )
334  {
335  *y = alpha + (*x);
336  x += incx;
337  y += incy;
338  }
339  }
340  }
341 
342  template LIB_UTILITIES_EXPORT void Sadd( int n, const Nektar::NekDouble alpha, const Nektar::NekDouble *x,
343  const int incx, Nektar::NekDouble *y, const int incy);
344 
345  /// \brief Subtract vector z = x-y
346  template<class T> void Vsub( int n, const T *x, const int incx, const T *y,
347  const int incy, T *z, const int incz)
348  {
349  ++n;
350  if (incx == 1 && incy == 1 && incz == 1)
351  {
352  while( --n )
353  {
354  *z = (*x) - (*y);
355  ++x;
356  ++y;
357  ++z;
358  }
359  }
360  else
361  {
362  while( --n )
363  {
364  *z = (*x) - (*y);
365  x += incx;
366  y += incy;
367  z += incz;
368  }
369  }
370  }
371 
372  template LIB_UTILITIES_EXPORT void Vsub( int n, const Nektar::NekDouble *x, const int incx, const Nektar::NekDouble *y,
373  const int incy, Nektar::NekDouble *z, const int incz);
374 
375  /// \brief Zero vector
376  template<class T> void Zero(int n, T *x, const int incx)
377  {
378  if(incx == 1)
379  {
380  std::memset(x,'\0', n*sizeof(T));
381  }
382  else
383  {
384  T zero = 0;
385  ++n;
386  while(--n)
387  {
388  *x = zero;
389  x+=incx;
390  }
391  }
392  }
393 
394  template LIB_UTILITIES_EXPORT void Zero(int n, Nektar::NekDouble *x, const int incx);
395  template LIB_UTILITIES_EXPORT void Zero(int n, int *x, const int incx);
396  template LIB_UTILITIES_EXPORT void Zero(int n, long *x, const int incx);
397 
398  /// \brief Negate x = -x
399  template<class T> void Neg( int n, T *x, const int incx)
400  {
401  while( n-- )
402  {
403  *x = -(*x);
404  x += incx;
405  }
406  }
407 
408  template LIB_UTILITIES_EXPORT void Neg( int n, Nektar::NekDouble *x, const int incx);
409 
410  /// \brief sqrt y = sqrt(x)
411  template<class T> void Vsqrt(int n, const T *x, const int incx,
412  T *y, const int incy)
413  {
414  while (n--)
415  {
416  *y = sqrt( *x );
417  x += incx;
418  y += incy;
419  }
420  }
421 
422  template LIB_UTILITIES_EXPORT void Vsqrt(int n, const Nektar::NekDouble *x, const int incx,
423  Nektar::NekDouble *y, const int incy);
424 
425 
426  /// \brief vabs: y = |x|
427  template<class T> void Vabs(int n, const T *x, const int incx,
428  T *y, const int incy)
429  {
430  while( n-- )
431  {
432  *y = ( *x >0)? *x:-(*x);
433  x += incx;
434  y += incy;
435  }
436  }
437 
438  template LIB_UTILITIES_EXPORT void Vabs(int n, const Nektar::NekDouble *x, const int incx,
439  Nektar::NekDouble *y, const int incy);
440 
441 
442  /********** Triad routines ***********************/
443 
444  /// \brief vvtvp (vector times vector plus vector): z = w*x + y
445  template<class T> void Vvtvp(int n,
446  const T *w, const int incw,
447  const T *x, const int incx,
448  const T *y, const int incy,
449  T *z, const int incz)
450  {
451  while( n-- )
452  {
453  *z = (*w) * (*x) + (*y);
454  w += incw;
455  x += incx;
456  y += incy;
457  z += incz;
458  }
459  }
460 
461  template LIB_UTILITIES_EXPORT void Vvtvp(int n,
462  const Nektar::NekDouble *w, const int incw,
463  const Nektar::NekDouble *x, const int incx,
464  const Nektar::NekDouble *y, const int incy,
465  Nektar::NekDouble *z, const int incz);
466 
467  /// \brief vvtvm (vector times vector plus vector): z = w*x - y
468  template<class T> void Vvtvm(int n, const T *w, const int incw, const T *x,
469  const int incx, const T *y, const int incy,
470  T *z, const int incz)
471  {
472  while( n-- )
473  {
474  *z = (*w) * (*x) - (*y);
475  w += incw;
476  x += incx;
477  y += incy;
478  z += incz;
479  }
480  }
481 
482  template LIB_UTILITIES_EXPORT void Vvtvm(int n, const Nektar::NekDouble *w, const int incw, const Nektar::NekDouble *x,
483  const int incx, const Nektar::NekDouble *y, const int incy,
484  Nektar::NekDouble *z, const int incz);
485 
486 
487  /// \brief svtvp (scalar times vector plus vector): z = alpha*x + y
488  template<class T> LIB_UTILITIES_EXPORT void Svtvp(int n, const T alpha, const T *x,
489  const int incx, const T *y, const int incy,
490  T *z, const int incz)
491  {
492  ++n;
493  if (incx == 1 && incy == 1 && incz == 1)
494  {
495  while( --n )
496  {
497  *z = alpha * (*x) + (*y);
498  ++x;
499  ++y;
500  ++z;
501  }
502  }
503  else
504  {
505  while( --n )
506  {
507  *z = alpha * (*x) + (*y);
508  x += incx;
509  y += incy;
510  z += incz;
511  }
512  }
513  }
514 
515  template LIB_UTILITIES_EXPORT void Svtvp(int n, const Nektar::NekDouble alpha, const Nektar::NekDouble *x,
516  const int incx, const Nektar::NekDouble *y, const int incy,
517  Nektar::NekDouble *z, const int incz);
518 
519 
520  /// \brief svtvp (scalar times vector plus vector): z = alpha*x - y
521  template<class T> void Svtvm(int n, const T alpha, const T *x,
522  const int incx, const T *y, const int incy,
523  T *z, const int incz)
524  {
525  while( n-- )
526  {
527  *z = alpha * (*x) - (*y);
528  x += incx;
529  y += incy;
530  z += incz;
531  }
532  }
533 
534  template LIB_UTILITIES_EXPORT void Svtvm(int n, const Nektar::NekDouble alpha, const Nektar::NekDouble *x,
535  const int incx, const Nektar::NekDouble *y, const int incy,
536  Nektar::NekDouble *z, const int incz);
537 
538  /// \brief vvtvvtp (vector times vector plus vector times vector):
539  // z = v*w + x*y
540  template<class T> void Vvtvvtp (int n,
541  const T* v, int incv,
542  const T* w, int incw,
543  const T* x, int incx,
544  const T* y, int incy,
545  T* z, int incz)
546  {
547  while( n-- )
548  {
549  *z = (*v) * (*w) + (*x) * (*y);
550  v += incv;
551  w += incw;
552  x += incx;
553  y += incy;
554  z += incz;
555  }
556  }
557  template LIB_UTILITIES_EXPORT void Vvtvvtp (int n,
558  const Nektar::NekDouble* v, int incv,
559  const Nektar::NekDouble* w, int incw,
560  const Nektar::NekDouble* x, int incx,
561  const Nektar::NekDouble* y, int incy,
562  Nektar::NekDouble* z, int incz);
563 
564 
565  /// \brief vvtvvtm (vector times vector minus vector times vector):
566  // z = v*w - x*y
567  template<class T> void Vvtvvtm (int n,
568  const T* v, int incv,
569  const T* w, int incw,
570  const T* x, int incx,
571  const T* y, int incy,
572  T* z, int incz)
573  {
574  while( n-- )
575  {
576  *z = (*v) * (*w) - (*x) * (*y);
577  v += incv;
578  w += incw;
579  x += incx;
580  y += incy;
581  z += incz;
582  }
583  }
584 
585  template LIB_UTILITIES_EXPORT void Vvtvvtm (int n,
586  const Nektar::NekDouble* v, int incv,
587  const Nektar::NekDouble* w, int incw,
588  const Nektar::NekDouble* x, int incx,
589  const Nektar::NekDouble* y, int incy,
590  Nektar::NekDouble* z, int incz);
591 
592  /// \brief vvtvvtp (scalar times vector plus scalar times vector):
593  // z = alpha*x + beta*y
594  template<class T> void Svtsvtp (int n,
595  const T alpha,
596  const T* x, int incx,
597  const T beta,
598  const T* y, int incy,
599  T* z, int incz)
600  {
601  while( n-- )
602  {
603  *z = alpha * (*x) + beta * (*y);
604  x += incx;
605  y += incy;
606  z += incz;
607  }
608  }
609 
610  template LIB_UTILITIES_EXPORT void Svtsvtp (int n,
611  const Nektar::NekDouble alpha,
612  const Nektar::NekDouble* x, int incx,
613  const Nektar::NekDouble beta,
614  const Nektar::NekDouble* y, int incy,
615  Nektar::NekDouble* z, int incz);
616 
617 
618  /// \brief Vstvpp (scalar times vector plus vector plus vector):
619  // z = v*w + x*y
620  template<class T> void Vstvpp(int n,
621  const T alpha,
622  const T* v, int incv,
623  const T* w, int incw,
624  const T* x, int incx,
625  T* z, int incz)
626  {
627  while( n-- )
628  {
629  *z = alpha * (*v) + (*w) + (*x);
630  v += incv;
631  w += incw;
632  x += incx;
633  z += incz;
634  }
635  }
636 
637  template LIB_UTILITIES_EXPORT void Vstvpp(int n,
638  const Nektar::NekDouble alpha,
639  const Nektar::NekDouble* v, int incv,
640  const Nektar::NekDouble* w, int incw,
641  const Nektar::NekDouble* x, int incx,
642  Nektar::NekDouble* z, int incz);
643 
644  /************ Misc routine from Veclib (and extras) ************/
645 
646  /// \brief Gather vector z[i] = x[y[i]]
647  template<class T> void Gathr(int n, const T *x, const int *y,
648  T *z)
649  {
650  while (n--)
651  {
652  *z++ = *(x + *y++);
653  }
654  return;
655  }
656 
657  template LIB_UTILITIES_EXPORT void Gathr(int n, const Nektar::NekDouble *x, const int *y,
658  Nektar::NekDouble *z);
659 
660 
661  /// \brief Gather vector z[i] = sign[i]*x[y[i]]
662  template<class T> void Gathr(int n, const T *sign, const T *x, const int *y,
663  T *z)
664  {
665  while (n--)
666  {
667  *z++ = *(sign++) * (*(x + *y++));
668  }
669  return;
670  }
671 
672  template LIB_UTILITIES_EXPORT void Gathr(int n, const Nektar::NekDouble *sign, const Nektar::NekDouble *x, const int *y,
673  Nektar::NekDouble *z);
674 
675  /// \brief Scatter vector z[y[i]] = x[i]
676  template<class T> void Scatr(int n, const T *x, const int *y,
677  T *z)
678  {
679  while (n--)
680  {
681  *(z + *(y++)) = *(x++);
682  }
683  }
684 
685  template LIB_UTILITIES_EXPORT void Scatr(int n, const Nektar::NekDouble *x, const int *y,
686  Nektar::NekDouble *z);
687 
688  /// \brief Scatter vector z[y[i]] = sign[i]*x[i]
689  template<class T> void Scatr(int n, const T *sign, const T *x, const int *y,
690  T *z)
691  {
692  while (n--)
693  {
694  if(*sign)
695  {
696  *(z + *(y++)) = *(sign++) * (*(x++));
697  }
698  else
699  {
700  x++;
701  y++;
702  sign++;
703  }
704  }
705  }
706 
707  template LIB_UTILITIES_EXPORT void Scatr(int n, const Nektar::NekDouble *sign, const Nektar::NekDouble *x, const int *y,
708  Nektar::NekDouble *z);
709 
710 
711  /// \brief Assemble z[y[i]] += x[i]; z should be zero'd first
712  template<class T> void Assmb(int n, const T *x, const int *y,
713  T *z)
714  {
715  while (n--)
716  {
717  *(z + *(y++)) += *(x++);
718  }
719  }
720 
721  template LIB_UTILITIES_EXPORT void Assmb(int n, const Nektar::NekDouble *x, const int *y,
722  Nektar::NekDouble *z);
723 
724  /// \brief Assemble z[y[i]] += sign[i]*x[i]; z should be zero'd first
725  template<class T> void Assmb(int n, const T *sign, const T *x, const int *y,
726  T *z)
727  {
728  while (n--)
729  {
730  *(z + *(y++)) += *(sign++) * (*(x++));
731  }
732  }
733 
734  template LIB_UTILITIES_EXPORT void Assmb(int n, const Nektar::NekDouble *sign, const Nektar::NekDouble *x, const int *y,
735  Nektar::NekDouble *z);
736 
737  /************* Reduction routines *****************/
738 
739  /// \brief Subtract return sum(x)
740  template<class T> T Vsum( int n, const T *x, const int incx)
741  {
742 
743  T sum = 0;
744 
745  while( n-- )
746  {
747  sum += (*x);
748  x += incx;
749  }
750 
751  return sum;
752  }
753 
754  template LIB_UTILITIES_EXPORT Nektar::NekDouble Vsum( int n, const Nektar::NekDouble *x, const int incx);
755  template LIB_UTILITIES_EXPORT int Vsum( int n, const int *x, const int incx);
756 
757  /// \brief Return the index of the maximum element in x
758  template<class T> int Imax( int n, const T *x, const int incx)
759  {
760 
761  int i, indx = ( n > 0 ) ? 0 : -1;
762  T xmax = *x;
763 
764  for (i = 0; i < n; i++)
765  {
766  if (*x > xmax)
767  {
768  xmax = *x;
769  indx = i;
770  }
771  x += incx;
772  }
773 
774  return indx;
775  }
776 
777  template LIB_UTILITIES_EXPORT int Imax( int n, const Nektar::NekDouble *x, const int incx);
778  template LIB_UTILITIES_EXPORT int Imax( int n, const int *x, const int incx);
779 
780  /// \brief Return the maximum element in x -- called vmax to avoid
781  /// conflict with max
782  template<class T> T Vmax( int n, const T *x, const int incx)
783  {
784 
785  T xmax = *x;
786 
787  while( n-- )
788  {
789  if (*x > xmax)
790  {
791  xmax = *x;
792  }
793  x += incx;
794  }
795 
796  return xmax;
797  }
798 
799  template LIB_UTILITIES_EXPORT Nektar::NekDouble Vmax( int n, const Nektar::NekDouble *x, const int incx);
800  template LIB_UTILITIES_EXPORT int Vmax( int n, const int *x, const int incx);
801 
802  /// \brief Return the index of the maximum absolute element in x
803  template<class T> int Iamax( int n, const T *x, const int incx)
804  {
805 
806  int i, indx = ( n > 0 ) ? 0 : -1;
807  T xmax = *x;
808  T xm;
809 
810  for (i = 0; i < n; i++)
811  {
812  xm = (*x > 0)? *x: -*x;
813  if (xm > xmax)
814  {
815  xmax = xm;
816  indx = i;
817  }
818  x += incx;
819  }
820 
821  return indx;
822  }
823 
824  template LIB_UTILITIES_EXPORT int Iamax( int n, const Nektar::NekDouble *x, const int incx);
825 
826  /// \brief Return the maximum absolute element in x
827  /// called vamax to avoid conflict with max
828  template<class T> T Vamax( int n, const T *x, const int incx)
829  {
830 
831  T xmax = *x;
832  T xm;
833 
834  while( n-- )
835  {
836  xm = (*x > 0)? *x: -*x;
837  if (xm > xmax)
838  {
839  xmax = xm;
840  }
841  x += incx;
842  }
843  return xmax;
844  }
845 
846  template LIB_UTILITIES_EXPORT Nektar::NekDouble Vamax( int n, const Nektar::NekDouble *x, const int incx);
847 
848 
849  /// \brief Return the index of the minimum element in x
850  template<class T> int Imin( int n, const T *x, const int incx)
851  {
852 
853  int i, indx = ( n > 0 ) ? 0 : -1;
854  T xmin = *x;
855 
856  for(i = 0;i < n;i++)
857  {
858  if( *x < xmin )
859  {
860  xmin = *x;
861  indx = i;
862  }
863  x += incx;
864  }
865 
866  return indx;
867  }
868 
869  template LIB_UTILITIES_EXPORT int Imin( int n, const Nektar::NekDouble *x, const int incx);
870  template LIB_UTILITIES_EXPORT int Imin( int n, const int *x, const int incx);
871 
872  /// \brief Return the minimum element in x - called vmin to avoid
873  /// conflict with min
874  template<class T> T Vmin( int n, const T *x, const int incx)
875  {
876 
877  T xmin = *x;
878 
879  while( n-- )
880  {
881  if (*x < xmin)
882  {
883  xmin = *x;
884  }
885  x += incx;
886  }
887 
888  return xmin;
889  }
890 
891  template LIB_UTILITIES_EXPORT Nektar::NekDouble Vmin( int n, const Nektar::NekDouble *x, const int incx);
892  template LIB_UTILITIES_EXPORT int Vmin( int n, const int *x, const int incx);
893 
894  /// \brief Return number of NaN elements of x
895  template<class T> int Nnan(int n, const T *x, const int incx)
896  {
897 
898  int nNan = 0;
899 
900  while (n--)
901  {
902  if (*x != *x)
903  {
904  nNan++;
905  }
906  x += incx;
907  }
908 
909  return nNan;
910  }
911 
912  template LIB_UTILITIES_EXPORT int Nnan(int n, const Nektar::NekDouble *x, const int incx);
913  template LIB_UTILITIES_EXPORT int Nnan(int n, const float *x, const int incx);
914  template LIB_UTILITIES_EXPORT int Nnan(int n, const int *x, const int incx);
915 
916  /// \brief vvtvp (vector times vector times vector): z = w*x*y
917  template<class T> T Dot( int n,
918  const T *w,
919  const T *x)
920  {
921  T sum = 0;
922 
923  while( n-- )
924  {
925  sum += (*w) * (*x);
926  ++w;
927  ++x;
928  }
929  return sum;
930  }
931 
932  template LIB_UTILITIES_EXPORT Nektar::NekDouble Dot( int n,
933  const Nektar::NekDouble *w,
934  const Nektar::NekDouble *x);
935 
936  /// \brief vvtvp (vector times vector times vector): z = w*x*y
937  template<class T> T Dot( int n,
938  const T *w, const int incw,
939  const T *x, const int incx)
940  {
941  T sum = 0;
942 
943  while( n-- )
944  {
945  sum += (*w) * (*x);
946  w += incw;
947  x += incx;
948  }
949  return sum;
950  }
951 
952  template LIB_UTILITIES_EXPORT Nektar::NekDouble Dot( int n,
953  const Nektar::NekDouble *w, const int incw,
954  const Nektar::NekDouble *x, const int incx);
955 
956  /// \brief vvtvp (vector times vector times vector): z = w*x*y
957  template<class T> T Dot2( int n,
958  const T *w,
959  const T *x,
960  const int *y)
961  {
962  T sum = 0;
963 
964  while( n-- )
965  {
966  sum += (*y == 1 ? (*w) * (*x) : 0 );
967  ++w;
968  ++x;
969  ++y;
970  }
971  return sum;
972  }
973 
974  template LIB_UTILITIES_EXPORT Nektar::NekDouble Dot2( int n,
975  const Nektar::NekDouble *w,
976  const Nektar::NekDouble *x,
977  const int *y);
978 
979 
980  /// \brief vvtvp (vector times vector times vector): z = w*x*y
981  template<class T> T Dot2( int n,
982  const T *w, const int incw,
983  const T *x, const int incx,
984  const int *y, const int incy)
985  {
986  T sum = 0;
987 
988  while( n-- )
989  {
990  sum += (*y == 1 ? (*w) * (*x) : 0.0 );
991  w += incw;
992  x += incx;
993  y += incy;
994  }
995  return sum;
996  }
997 
998  template LIB_UTILITIES_EXPORT Nektar::NekDouble Dot2( int n,
999  const Nektar::NekDouble *w, const int incw,
1000  const Nektar::NekDouble *x, const int incx,
1001  const int *y, const int incy);
1002 
1003 /*
1004  // \brief copy one int vector to another
1005  void Vcopy(int n, const int *x, const int incx, int *y,
1006  const int incy)
1007  {
1008  if( incx ==1 && incy == 1)
1009  {
1010  memcpy(y,x,n*sizeof(int));
1011  }
1012  else
1013  {
1014  while( n-- )
1015  {
1016  *y = *x;
1017  x += incx;
1018  y += incy;
1019  }
1020  }
1021  }
1022 
1023  // \brief copy one unsigned int vector to another
1024  void Vcopy(int n, const unsigned int *x, const int incx, unsigned int *y,
1025  const int incy)
1026  {
1027  if( incx ==1 && incy == 1)
1028  {
1029  memcpy(y,x,n*sizeof(int));
1030  }
1031  else
1032  {
1033  while( n-- )
1034  {
1035  *y = *x;
1036  x += incx;
1037  y += incy;
1038  }
1039  }
1040  }
1041 
1042  // \brief copy one double vector to another
1043  void Vcopy(int n, const double *x, const int incx, double *y,
1044  const int incy)
1045  {
1046  if( incx ==1 && incy == 1)
1047  {
1048  memcpy(y,x,n*sizeof(double));
1049  }
1050  else
1051  {
1052  while( n-- )
1053  {
1054  *y = *x;
1055  x += incx;
1056  y += incy;
1057  }
1058  }
1059  }
1060 */
1061 
1062  // \brief copy one vector to another
1063  template<typename T>
1064  void Vcopy(int n, const T *x, const int incx,
1065  T *y, const int incy)
1066  {
1067  if( incx ==1 && incy == 1)
1068  {
1069  memcpy(y,x,n*sizeof(T));
1070  }
1071  else
1072  {
1073  while( n-- )
1074  {
1075  *y = *x;
1076  x += incx;
1077  y += incy;
1078  }
1079  }
1080  }
1081 
1082  template LIB_UTILITIES_EXPORT void Vcopy( int n, const int *x, const int incx, int *y, const int incy);
1083  template LIB_UTILITIES_EXPORT void Vcopy( int n, const unsigned int *x, const int incx, unsigned int *y, const int incy);
1084  template LIB_UTILITIES_EXPORT void Vcopy( int n, const Nektar::NekDouble *x, const int incx, Nektar::NekDouble *y, const int incy);
1085 
1086 
1087  // \brief reverse the ordering of vector to another
1088  template<class T> void Reverse( int n, const T *x, const int incx, T *y, const int incy)
1089  {
1090  int i;
1091  T store;
1092 
1093  // Perform element by element swaps in case x and y reference the same
1094  // array.
1095  int nloop = n/2;
1096 
1097  // copy value in case of n is odd number
1098  y[nloop] = x[nloop];
1099 
1100  const T* x_end = x + (n-1)*incx;
1101  T* y_end = y + (n-1)*incy;
1102  for (i = 0; i < nloop; ++i) {
1103  store = *x_end;
1104  *y_end = *x;
1105  *y = store;
1106  x += incx;
1107  y += incy;
1108  x_end -= incx;
1109  y_end -= incy;
1110  }
1111  }
1112 
1113  template LIB_UTILITIES_EXPORT void Reverse( int n, const Nektar::NekDouble *x, const int incx, Nektar::NekDouble *y, const int incy);
1114 }
IMM1
#define IMM1
Definition: Vmath.cpp:59
Vmath
Definition: Vmath.cpp:39
IM1
#define IM1
Definition: Vmath.cpp:56
IM2
#define IM2
Definition: Vmath.cpp:57
Vmath::Gathr
void Gathr(int n, const T *x, const int *y, T *z)
Gather vector z[i] = x[y[i]].
Definition: Vmath.cpp:647
sign
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:16
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:411
AM
#define AM
Definition: Vmath.cpp:58
Vmath::Vmax
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:782
Vmath::Vmin
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
Definition: Vmath.cpp:874
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
Vmath::Imin
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
Definition: Vmath.cpp:850
Vmath::Svtvp
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:488
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:445
Vmath::Sdiv
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/y.
Definition: Vmath.cpp:274
IQ1
#define IQ1
Definition: Vmath.cpp:62
Vmath::Vdiv
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:244
Vmath::Vabs
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
Definition: Vmath.cpp:427
Vmath::Imax
int Imax(int n, const T *x, const int incx)
Return the index of the maximum element in x.
Definition: Vmath.cpp:758
NektarUnivTypeDefs.hpp
IR1
#define IR1
Definition: Vmath.cpp:64
Vmath::Reverse
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1088
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:216
IQ2
#define IQ2
Definition: Vmath.cpp:63
NTAB
#define NTAB
Definition: Vmath.cpp:66
Vmath::ran2
T ran2(long *idum)
Generates a number from ~Normal(0,1)
Definition: Vmath.cpp:72
LIB_UTILITIES_EXPORT
#define LIB_UTILITIES_EXPORT
Definition: LibUtilitiesDeclspec.h:42
IR2
#define IR2
Definition: Vmath.cpp:65
Vmath::Nnan
int Nnan(int n, const T *x, const int incx)
Return number of NaN elements of x.
Definition: Vmath.cpp:895
IA2
#define IA2
Definition: Vmath.cpp:61
Vmath::Svtvm
void Svtvm(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x - y
Definition: Vmath.cpp:521
Vmath::Scatr
void Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
Definition: Vmath.cpp:676
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:399
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Assmb
void Assmb(int n, const T *x, const int *y, T *z)
Assemble z[y[i]] += x[i]; z should be zero&#39;d first.
Definition: Vmath.cpp:712
Vmath::Vvtvvtm
void Vvtvvtm(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtm (vector times vector minus vector times vector):
Definition: Vmath.cpp:567
Vmath::Vamax
T Vamax(int n, const T *x, const int incx)
Return the maximum absolute element in x called vamax to avoid conflict with max. ...
Definition: Vmath.cpp:828
Vmath::Iamax
int Iamax(int n, const T *x, const int incx)
Return the index of the maximum absolute element in x.
Definition: Vmath.cpp:803
IA1
#define IA1
Definition: Vmath.cpp:60
Vmath::Sadd
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.cpp:318
LibUtilitiesDeclspec.h
Vmath::Dot
T Dot(int n, const T *w, const T *x)
vvtvp (vector times vector times vector): z = w*x*y
Definition: Vmath.cpp:917
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:346
RNMX
#define RNMX
Definition: Vmath.cpp:69
Vmath::Vstvpp
void Vstvpp(int n, const T alpha, const T *v, int incv, const T *w, int incw, const T *x, int incx, T *z, int incz)
Vstvpp (scalar times vector plus vector plus vector):
Definition: Vmath.cpp:620
Vmath.hpp
Vmath::Vvtvvtp
void Vvtvvtp(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtp (vector times vector plus vector times vector):
Definition: Vmath.cpp:540
Vmath::Vvtvm
void Vvtvm(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvm (vector times vector plus vector): z = w*x - y
Definition: Vmath.cpp:468
Vmath::Dot2
T Dot2(int n, const T *w, const T *x, const int *y)
vvtvp (vector times vector times vector): z = w*x*y
Definition: Vmath.cpp:957
Vmath::FillWhiteNoise
void FillWhiteNoise(int n, const T eps, T *x, const int incx, int outseed)
Fills a vector with white noise.
Definition: Vmath.cpp:139
Vmath::Svtsvtp
void Svtsvtp(int n, const T alpha, const T *x, int incx, const T beta, const T *y, int incy, T *z, int incz)
vvtvvtp (scalar times vector plus scalar times vector):
Definition: Vmath.cpp:594
NDIV
#define NDIV
Definition: Vmath.cpp:67
Vmath::Vsum
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:740
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:376
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:302
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:186
Generated on Wed Jan 20 2021 03:33:43 for Nektar++ by   doxygen 1.8.13