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 
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  template LIB_UTILITIES_EXPORT void Fill( int n,
56  const Nektar::NekSingle alpha, Nektar::NekSingle *x, const int incx );
57 
58  #define IM1 2147483563
59  #define IM2 2147483399
60  #define AM (1.0/IM1)
61  #define IMM1 (IM1-1)
62  #define IA1 40014
63  #define IA2 40692
64  #define IQ1 53668
65  #define IQ2 52774
66  #define IR1 12211
67  #define IR2 3791
68  #define NTAB 32
69  #define NDIV (1+IMM1/NTAB)
70  #define EPS 1.2e-7
71  #define RNMX (1.0-EPS)
72 
73  /// \brief Generates a number from ~Normal(0,1)
74  template<class T> T ran2 (long* idum)
75  /* ------------------------------------------------------------------------- *
76  * Ran2 from NR 2e. Returns a uniform random deviate between 0.0 &
77  * 1.0 (exclusive of endpoints). Call with idum a negative integer to
78  * initialize; thereafter, do not alter idum between successive
79  * deviates in a sequence. RNMX should approximate the largest
80  * floating value that is less than 1.
81  * ------------------------------------------------------------------------- */
82  {
83  int j;
84  long k;
85  static long idum2=123456789;
86  static long iy=0;
87  static long iv[NTAB];
88  T temp;
89 
90  if (*idum <= 0) {
91  if (-(*idum) < 1) *idum = 1;
92  else *idum = -(*idum);
93  idum2 = (*idum);
94  for (j=NTAB+7; j>=0; j--) {
95  k = (*idum) / IQ1;
96  *idum = IA1 * (*idum - k*IQ1) - k*IR1;
97  if (*idum < 0) *idum += IM1;
98  if (j < NTAB) iv[j] = *idum;
99  }
100  iy = iv[0];
101  }
102 
103  k = (*idum) / IQ1;
104  *idum = IA1*(*idum - k*IQ1) - k*IR1;
105  if (*idum < 0) *idum += IM1;
106 
107  k = idum2 / IQ2;
108  idum2 = IA2*(idum2 - k*IQ2) - k*IR2;
109  if (idum2 < 0) idum2 += IM2;
110 
111  j = iy / NDIV;
112  iy = iv[j] - idum2;
113  iv[j] = *idum;
114  if (iy < 1) iy += IMM1;
115 
116  if ((temp=AM*iy) > RNMX) return RNMX;
117  else return temp;
118  }
119 
120  #undef IM1
121  #undef IM2
122  #undef AM
123  #undef IMM1
124  #undef IA1
125  #undef IA2
126  #undef IQ1
127  #undef IQ2
128  #undef IR1
129  #undef IR2
130  #undef NTAB
131  #undef NDIV
132  #undef EPS
133  #undef RNMX
134 
135 #ifdef NEKTAR_USE_THREAD_SAFETY
136  static boost::mutex mutex;
137 #endif
138  template LIB_UTILITIES_EXPORT Nektar::NekDouble ran2 (long* idum);
139  template LIB_UTILITIES_EXPORT Nektar::NekSingle ran2 (long* idum);
140 
141  /// \brief Fills a vector with white noise.
142  template<class T> void FillWhiteNoise( int n, const T eps, T* x,
143  const int incx, int outseed)
144  {
145 #ifdef NEKTAR_USE_THREAD_SAFETY
146  // Protect the static vars here and in ran2
147  boost::mutex::scoped_lock l(mutex);
148 #endif
149 
150  // Define static variables for generating random numbers
151  static int iset = 0;
152  static T gset;
153  static long seed = 0;
154 
155  // Bypass seed if outseed was specified
156  if( outseed != 9999)
157  {
158  seed = long(outseed);
159  }
160 
161  while( n-- )
162  {
163  T fac, rsq, v1, v2;
164 
165  if (iset == 0)
166  {
167  do
168  {
169  v1 = 2.0 * ran2<T> (&seed) - 1.0;
170  v2 = 2.0 * ran2<T> (&seed) - 1.0;
171  rsq = v1*v1 + v2*v2;
172  } while (rsq >= 1.0 || rsq == 0.0);
173  fac = sqrt(-2.0 * log (rsq) / rsq);
174  gset = v1 * fac;
175  iset = 1;
176  *x = eps * v2 * fac;
177  }
178  else
179  {
180  iset = 0;
181  *x = eps * gset;
182  }
183  x += incx;
184  }
185  }
186  template LIB_UTILITIES_EXPORT void FillWhiteNoise( int n, const Nektar::NekDouble eps, Nektar::NekDouble *x, const int incx, int outseed);
187  template LIB_UTILITIES_EXPORT void FillWhiteNoise( int n,
188  const Nektar::NekSingle eps, Nektar::NekSingle *x,
189  const int incx, int outseed);
190 
191  /// \brief Multiply vector z = x*y
192  template<class T> void Vmul( int n, const T *x, const int incx, const T *y,
193  const int incy, T *z, const int incz)
194  {
195  ++n;
196  if (incx == 1 && incy == 1 && incz == 1)
197  {
198  while( --n )
199  {
200  *z = (*x) * (*y);
201  ++x;
202  ++y;
203  ++z;
204  }
205  }
206  else
207  {
208  while( --n )
209  {
210  *z = (*x) * (*y);
211  x += incx;
212  y += incy;
213  z += incz;
214  }
215  }
216  }
217  template LIB_UTILITIES_EXPORT void Vmul( int n, const Nektar::NekDouble *x, const int incx, const Nektar::NekDouble *y,
218  const int incy, Nektar::NekDouble*z, const int incz);
219  template LIB_UTILITIES_EXPORT void Vmul( int n, const Nektar::NekSingle *x,
220  const int incx, const Nektar::NekSingle *y,
221  const int incy, Nektar::NekSingle *z, const int incz);
222 
223  /// \brief Scalar multiply y = alpha*x
224 
225  template<class T> void Smul( int n, const T alpha, const T *x, const int incx,
226  T *y, const int incy)
227  {
228  ++n;
229  if (incx == 1 && incy == 1)
230  {
231  while( --n )
232  {
233  *y = alpha * (*x);
234  ++x;
235  ++y;
236  }
237  }
238  else
239  {
240  while( --n )
241  {
242  *y = alpha * (*x);
243  x += incx;
244  y += incy;
245  }
246  }
247  }
248 
249  template LIB_UTILITIES_EXPORT void Smul( int n, const Nektar::NekDouble alpha, const Nektar::NekDouble *x, const int incx,
250  Nektar::NekDouble *y, const int incy);
251 
252  template LIB_UTILITIES_EXPORT void Smul( int n,
253  const Nektar::NekSingle alpha, const Nektar::NekSingle *x,
254  const int incx, Nektar::NekSingle *y, const int incy);
255 
256  /// \brief Multiply vector z = x/y
257  template<class T> void Vdiv( int n, const T *x, const int incx, const T *y,
258  const int incy, T*z, const int incz)
259  {
260  ++n;
261  if (incx == 1 && incy == 1)
262  {
263  while( --n )
264  {
265  *z = (*x) / (*y);
266  ++x;
267  ++y;
268  ++z;
269  }
270  }
271  else
272  {
273  while( --n )
274  {
275  *z = (*x) / (*y);
276  x += incx;
277  y += incy;
278  z += incz;
279  }
280  }
281  }
282 
283  template LIB_UTILITIES_EXPORT void Vdiv( int n, const Nektar::NekDouble *x, const int incx, const Nektar::NekDouble *y,
284  const int incy, Nektar::NekDouble *z, const int incz);
285 
286  template LIB_UTILITIES_EXPORT void Vdiv( int n, const Nektar::NekSingle *x,
287  const int incx, const Nektar::NekSingle *y,
288  const int incy, Nektar::NekSingle *z, const int incz);
289 
290  /// \brief Scalar multiply y = alpha/y
291  template<class T> void Sdiv( int n, const T alpha, const T *x,
292  const int incx, T *y, const int incy)
293  {
294  ++n;
295  if (incx == 1 && incy == 1)
296  {
297  while( --n )
298  {
299  *y = alpha / (*x);
300  ++x;
301  ++y;
302  }
303  }
304  else
305  {
306  while( --n )
307  {
308  *y = alpha / (*x);
309  x += incx;
310  y += incy;
311  }
312  }
313  }
314 
315  template LIB_UTILITIES_EXPORT void Sdiv( int n, const Nektar::NekDouble alpha, const Nektar::NekDouble *x,
316  const int incx, Nektar::NekDouble *y, const int incy);
317  template LIB_UTILITIES_EXPORT void Sdiv( int n,
318  const Nektar::NekSingle alpha, const Nektar::NekSingle *x,
319  const int incx, Nektar::NekSingle *y, const int incy);
320 
321  /// \brief Add vector z = x+y
322  template<class T> void Vadd( int n, const T *x, const int incx, const T *y,
323  const int incy, T *z, const int incz)
324  {
325  while( n-- )
326  {
327  *z = (*x) + (*y);
328  x += incx;
329  y += incy;
330  z += incz;
331  }
332  }
333 
334  template LIB_UTILITIES_EXPORT void Vadd( int n, const Nektar::NekDouble *x, const int incx, const Nektar::NekDouble *y,
335  const int incy, Nektar::NekDouble *z, const int incz);
336  template LIB_UTILITIES_EXPORT void Vadd( int n,
337  const Nektar::NekSingle *x, const int incx, const Nektar::NekSingle *y,
338  const int incy, Nektar::NekSingle *z, const int incz);
339 
340  /// \brief Add vector y = alpha - x
341  template<class T> void Sadd( int n, const T alpha, const T *x,
342  const int incx, T *y, const int incy)
343  {
344  ++n;
345  if (incx == 1 && incy == 1)
346  {
347  while( --n )
348  {
349  *y = alpha + (*x);
350  ++x;
351  ++y;
352  }
353  }
354  else
355  {
356  while( --n )
357  {
358  *y = alpha + (*x);
359  x += incx;
360  y += incy;
361  }
362  }
363  }
364 
365  template LIB_UTILITIES_EXPORT void Sadd( int n, const Nektar::NekDouble alpha, const Nektar::NekDouble *x,
366  const int incx, Nektar::NekDouble *y, const int incy);
367  template LIB_UTILITIES_EXPORT void Sadd( int n,
368  const Nektar::NekSingle alpha, const Nektar::NekSingle *x,
369  const int incx, Nektar::NekSingle *y, const int incy);
370 
371  /// \brief Subtract vector z = x-y
372  template<class T> void Vsub( int n, const T *x, const int incx, const T *y,
373  const int incy, T *z, const int incz)
374  {
375  ++n;
376  if (incx == 1 && incy == 1 && incz == 1)
377  {
378  while( --n )
379  {
380  *z = (*x) - (*y);
381  ++x;
382  ++y;
383  ++z;
384  }
385  }
386  else
387  {
388  while( --n )
389  {
390  *z = (*x) - (*y);
391  x += incx;
392  y += incy;
393  z += incz;
394  }
395  }
396  }
397 
398  template LIB_UTILITIES_EXPORT void Vsub( int n, const Nektar::NekDouble *x, const int incx, const Nektar::NekDouble *y,
399  const int incy, Nektar::NekDouble *z, const int incz);
400  template LIB_UTILITIES_EXPORT void Vsub( int n,
401  const Nektar::NekSingle *x, const int incx, const Nektar::NekSingle *y,
402  const int incy, Nektar::NekSingle *z, const int incz);
403 
404  /// \brief Add vector y = alpha - x
405  template<class T> void Ssub( int n, const T alpha, const T *x,
406  const int incx, T *y, const int incy)
407  {
408  ++n;
409  if (incx == 1 && incy == 1)
410  {
411  while( --n )
412  {
413  *y = alpha - (*x);
414  ++x;
415  ++y;
416  }
417  }
418  else
419  {
420  while( --n )
421  {
422  *y = alpha - (*x);
423  x += incx;
424  y += incy;
425  }
426  }
427  }
428 
429  template LIB_UTILITIES_EXPORT void Ssub( int n, const Nektar::NekDouble alpha, const Nektar::NekDouble *x,
430  const int incx, Nektar::NekDouble *y, const int incy);
431  template LIB_UTILITIES_EXPORT void Ssub( int n,
432  const Nektar::NekSingle alpha, const Nektar::NekSingle *x,
433  const int incx, Nektar::NekSingle *y, const int incy);
434 
435  /// \brief Zero vector
436  template<class T> void Zero(int n, T *x, const int incx)
437  {
438  if(incx == 1)
439  {
440  std::memset(x,'\0', n*sizeof(T));
441  }
442  else
443  {
444  T zero = 0;
445  ++n;
446  while(--n)
447  {
448  *x = zero;
449  x+=incx;
450  }
451  }
452  }
453 
454  template LIB_UTILITIES_EXPORT void Zero(int n, Nektar::NekDouble *x, const int incx);
455  template LIB_UTILITIES_EXPORT void Zero(int n, Nektar::NekSingle *x,
456  const int incx);
457  template LIB_UTILITIES_EXPORT void Zero(int n, int *x, const int incx);
458  template LIB_UTILITIES_EXPORT void Zero(int n, long *x, const int incx);
459 
460  /// \brief Negate x = -x
461  template<class T> void Neg( int n, T *x, const int incx)
462  {
463  while( n-- )
464  {
465  *x = -(*x);
466  x += incx;
467  }
468  }
469 
470  template LIB_UTILITIES_EXPORT void Neg( int n, Nektar::NekDouble *x, const int incx);
471  template LIB_UTILITIES_EXPORT void Neg( int n, Nektar::NekSingle *x,
472  const int incx);
473 
474  /// \brief sqrt y = sqrt(x)
475  template<class T> void Vsqrt(int n, const T *x, const int incx,
476  T *y, const int incy)
477  {
478  while (n--)
479  {
480  *y = sqrt( *x );
481  x += incx;
482  y += incy;
483  }
484  }
485 
486  template LIB_UTILITIES_EXPORT void Vsqrt(int n, const Nektar::NekDouble *x, const int incx,
487  Nektar::NekDouble *y, const int incy);
488  template LIB_UTILITIES_EXPORT void Vsqrt(int n, const Nektar::NekSingle *x,
489  const int incx, Nektar::NekSingle *y, const int incy);
490 
491 
492  /// \brief vabs: y = |x|
493  template<class T> void Vabs(int n, const T *x, const int incx,
494  T *y, const int incy)
495  {
496  while( n-- )
497  {
498  *y = ( *x >0)? *x:-(*x);
499  x += incx;
500  y += incy;
501  }
502  }
503 
504  template LIB_UTILITIES_EXPORT void Vabs(int n, const Nektar::NekDouble *x, const int incx,
505  Nektar::NekDouble *y, const int incy);
506  template LIB_UTILITIES_EXPORT void Vabs(int n, const Nektar::NekSingle *x,
507  const int incx, Nektar::NekSingle *y, const int incy);
508 
509 
510  /********** Triad routines ***********************/
511 
512  /// \brief vvtvp (vector times vector plus vector): z = w*x + y
513  template<class T> void Vvtvp(int n,
514  const T *w, const int incw,
515  const T *x, const int incx,
516  const T *y, const int incy,
517  T *z, const int incz)
518  {
519  while( n-- )
520  {
521  *z = (*w) * (*x) + (*y);
522  w += incw;
523  x += incx;
524  y += incy;
525  z += incz;
526  }
527  }
528 
529  template LIB_UTILITIES_EXPORT void Vvtvp(int n,
530  const Nektar::NekDouble *w, const int incw,
531  const Nektar::NekDouble *x, const int incx,
532  const Nektar::NekDouble *y, const int incy,
533  Nektar::NekDouble *z, const int incz);
534  template LIB_UTILITIES_EXPORT void Vvtvp(int n,
535  const Nektar::NekSingle *w, const int incw,
536  const Nektar::NekSingle *x, const int incx,
537  const Nektar::NekSingle *y, const int incy,
538  Nektar::NekSingle *z, const int incz);
539 
540  /// \brief vvtvm (vector times vector plus vector): z = w*x - y
541  template<class T> void Vvtvm(int n, const T *w, const int incw, const T *x,
542  const int incx, const T *y, const int incy,
543  T *z, const int incz)
544  {
545  while( n-- )
546  {
547  *z = (*w) * (*x) - (*y);
548  w += incw;
549  x += incx;
550  y += incy;
551  z += incz;
552  }
553  }
554 
555  template LIB_UTILITIES_EXPORT void Vvtvm(int n, const Nektar::NekDouble *w, const int incw, const Nektar::NekDouble *x,
556  const int incx, const Nektar::NekDouble *y, const int incy,
557  Nektar::NekDouble *z, const int incz);
558  template LIB_UTILITIES_EXPORT void Vvtvm(int n, const Nektar::NekSingle *w,
559  const int incw, const Nektar::NekSingle *x,
560  const int incx, const Nektar::NekSingle *y, const int incy,
561  Nektar::NekSingle *z, const int incz);
562 
563 
564  /// \brief svtvp (scalar times vector plus vector): z = alpha*x + y
565  template<class T> LIB_UTILITIES_EXPORT void Svtvp(int n, const T alpha, const T *x,
566  const int incx, const T *y, const int incy,
567  T *z, const int incz)
568  {
569  ++n;
570  if (incx == 1 && incy == 1 && incz == 1)
571  {
572  while( --n )
573  {
574  *z = alpha * (*x) + (*y);
575  ++x;
576  ++y;
577  ++z;
578  }
579  }
580  else
581  {
582  while( --n )
583  {
584  *z = alpha * (*x) + (*y);
585  x += incx;
586  y += incy;
587  z += incz;
588  }
589  }
590  }
591 
592  template LIB_UTILITIES_EXPORT void Svtvp(int n, const Nektar::NekDouble alpha, const Nektar::NekDouble *x,
593  const int incx, const Nektar::NekDouble *y, const int incy,
594  Nektar::NekDouble *z, const int incz);
595  template LIB_UTILITIES_EXPORT void Svtvp(int n,
596  const Nektar::NekSingle alpha, const Nektar::NekSingle *x,
597  const int incx, const Nektar::NekSingle *y, const int incy,
598  Nektar::NekSingle *z, const int incz);
599 
600 
601  /// \brief svtvp (scalar times vector plus vector): z = alpha*x - y
602  template<class T> void Svtvm(int n, const T alpha, const T *x,
603  const int incx, const T *y, const int incy,
604  T *z, const int incz)
605  {
606  while( n-- )
607  {
608  *z = alpha * (*x) - (*y);
609  x += incx;
610  y += incy;
611  z += incz;
612  }
613  }
614 
615  template LIB_UTILITIES_EXPORT void Svtvm(int n, const Nektar::NekDouble alpha, const Nektar::NekDouble *x,
616  const int incx, const Nektar::NekDouble *y, const int incy,
617  Nektar::NekDouble *z, const int incz);
618  template LIB_UTILITIES_EXPORT void Svtvm(int n,
619  const Nektar::NekSingle alpha, const Nektar::NekSingle *x,
620  const int incx, const Nektar::NekSingle *y, const int incy,
621  Nektar::NekSingle *z, const int incz);
622 
623  /// \brief vvtvvtp (vector times vector plus vector times vector):
624  // z = v*w + x*y
625  template<class T> void Vvtvvtp (int n,
626  const T* v, int incv,
627  const T* w, int incw,
628  const T* x, int incx,
629  const T* y, int incy,
630  T* z, int incz)
631  {
632  while( n-- )
633  {
634  *z = (*v) * (*w) + (*x) * (*y);
635  v += incv;
636  w += incw;
637  x += incx;
638  y += incy;
639  z += incz;
640  }
641  }
642  template LIB_UTILITIES_EXPORT void Vvtvvtp (int n,
643  const Nektar::NekDouble* v, int incv,
644  const Nektar::NekDouble* w, int incw,
645  const Nektar::NekDouble* x, int incx,
646  const Nektar::NekDouble* y, int incy,
647  Nektar::NekDouble* z, int incz);
648  template LIB_UTILITIES_EXPORT void Vvtvvtp (int n,
649  const Nektar::NekSingle* v, int incv,
650  const Nektar::NekSingle* w, int incw,
651  const Nektar::NekSingle* x, int incx,
652  const Nektar::NekSingle* y, int incy,
653  Nektar::NekSingle* z, int incz);
654 
655 
656  /// \brief vvtvvtm (vector times vector minus vector times vector):
657  // z = v*w - x*y
658  template<class T> void Vvtvvtm (int n,
659  const T* v, int incv,
660  const T* w, int incw,
661  const T* x, int incx,
662  const T* y, int incy,
663  T* z, int incz)
664  {
665  while( n-- )
666  {
667  *z = (*v) * (*w) - (*x) * (*y);
668  v += incv;
669  w += incw;
670  x += incx;
671  y += incy;
672  z += incz;
673  }
674  }
675 
676  template LIB_UTILITIES_EXPORT void Vvtvvtm (int n,
677  const Nektar::NekDouble* v, int incv,
678  const Nektar::NekDouble* w, int incw,
679  const Nektar::NekDouble* x, int incx,
680  const Nektar::NekDouble* y, int incy,
681  Nektar::NekDouble* z, int incz);
682  template LIB_UTILITIES_EXPORT void Vvtvvtm (int n,
683  const Nektar::NekSingle* v, int incv,
684  const Nektar::NekSingle* w, int incw,
685  const Nektar::NekSingle* x, int incx,
686  const Nektar::NekSingle* y, int incy,
687  Nektar::NekSingle* z, int incz);
688 
689  /// \brief vvtvvtp (scalar times vector plus scalar times vector):
690  // z = alpha*x + beta*y
691  template<class T> void Svtsvtp (int n,
692  const T alpha,
693  const T* x, int incx,
694  const T beta,
695  const T* y, int incy,
696  T* z, int incz)
697  {
698  while( n-- )
699  {
700  *z = alpha * (*x) + beta * (*y);
701  x += incx;
702  y += incy;
703  z += incz;
704  }
705  }
706 
707  template LIB_UTILITIES_EXPORT void Svtsvtp (int n,
708  const Nektar::NekDouble alpha,
709  const Nektar::NekDouble* x, int incx,
710  const Nektar::NekDouble beta,
711  const Nektar::NekDouble* y, int incy,
712  Nektar::NekDouble* z, int incz);
713  template LIB_UTILITIES_EXPORT void Svtsvtp (int n,
714  const Nektar::NekSingle alpha,
715  const Nektar::NekSingle* x, int incx,
716  const Nektar::NekSingle beta,
717  const Nektar::NekSingle* y, int incy,
718  Nektar::NekSingle* z, int incz);
719 
720 
721  /// \brief Vstvpp (scalar times vector plus vector plus vector):
722  // z = v*w + x*y
723  template<class T> void Vstvpp(int n,
724  const T alpha,
725  const T* v, int incv,
726  const T* w, int incw,
727  const T* x, int incx,
728  T* z, int incz)
729  {
730  while( n-- )
731  {
732  *z = alpha * (*v) + (*w) + (*x);
733  v += incv;
734  w += incw;
735  x += incx;
736  z += incz;
737  }
738  }
739 
740  template LIB_UTILITIES_EXPORT void Vstvpp(int n,
741  const Nektar::NekDouble alpha,
742  const Nektar::NekDouble* v, int incv,
743  const Nektar::NekDouble* w, int incw,
744  const Nektar::NekDouble* x, int incx,
745  Nektar::NekDouble* z, int incz);
746  template LIB_UTILITIES_EXPORT void Vstvpp(int n,
747  const Nektar::NekSingle alpha,
748  const Nektar::NekSingle* v, int incv,
749  const Nektar::NekSingle* w, int incw,
750  const Nektar::NekSingle* x, int incx,
751  Nektar::NekSingle* z, int incz);
752 
753  /************ Misc routine from Veclib (and extras) ************/
754 
755  /// \brief Gather vector z[i] = sign[i]*x[y[i]]
756  template<class T> void Gathr(int n, const T *sign, const T *x, const int *y,
757  T *z)
758  {
759  while (n--)
760  {
761  *z++ = *(sign++) * (*(x + *y++));
762  }
763  return;
764  }
765 
766  template LIB_UTILITIES_EXPORT void Gathr(int n, const Nektar::NekDouble *sign, const Nektar::NekDouble *x, const int *y,
767  Nektar::NekDouble *z);
768  template LIB_UTILITIES_EXPORT void Gathr(int n, const Nektar::NekSingle *sign, const Nektar::NekSingle *x, const int *y,
769  Nektar::NekSingle *z);
770 
771  /// \brief Scatter vector z[y[i]] = x[i]
772  template<class T> void Scatr(int n, const T *x, const int *y,
773  T *z)
774  {
775  while (n--)
776  {
777  *(z + *(y++)) = *(x++);
778  }
779  }
780 
781  template LIB_UTILITIES_EXPORT void Scatr(int n, const Nektar::NekDouble *x, const int *y,
782  Nektar::NekDouble *z);
783  template LIB_UTILITIES_EXPORT void Scatr(int n,
784  const Nektar::NekSingle *x, const int *y, Nektar::NekSingle *z);
785 
786  /// \brief Scatter vector z[y[i]] = sign[i]*x[i]
787  template<class T> void Scatr(int n, const T *sign, const T *x, const int *y,
788  T *z)
789  {
790  while (n--)
791  {
792  if(*sign)
793  {
794  *(z + *(y++)) = *(sign++) * (*(x++));
795  }
796  else
797  {
798  x++;
799  y++;
800  sign++;
801  }
802  }
803  }
804 
805  template LIB_UTILITIES_EXPORT void Scatr(int n, const Nektar::NekDouble *sign, const Nektar::NekDouble *x, const int *y,
806  Nektar::NekDouble *z);
807  template LIB_UTILITIES_EXPORT void Scatr(int n,
808  const Nektar::NekSingle *sign, const Nektar::NekSingle *x, const int *y,
809  Nektar::NekSingle *z);
810 
811 
812  /// \brief Assemble z[y[i]] += x[i]; z should be zero'd first
813  template<class T> void Assmb(int n, const T *x, const int *y,
814  T *z)
815  {
816  while (n--)
817  {
818  *(z + *(y++)) += *(x++);
819  }
820  }
821 
822  template LIB_UTILITIES_EXPORT void Assmb(int n, const Nektar::NekDouble *x, const int *y,
823  Nektar::NekDouble *z);
824  template LIB_UTILITIES_EXPORT void Assmb(int n, const Nektar::NekSingle *x,
825  const int *y, Nektar::NekSingle *z);
826 
827  /// \brief Assemble z[y[i]] += sign[i]*x[i]; z should be zero'd first
828  template<class T> void Assmb(int n, const T *sign, const T *x, const int *y,
829  T *z)
830  {
831  while (n--)
832  {
833  *(z + *(y++)) += *(sign++) * (*(x++));
834  }
835  }
836 
837  template LIB_UTILITIES_EXPORT void Assmb(int n, const Nektar::NekDouble *sign, const Nektar::NekDouble *x, const int *y,
838  Nektar::NekDouble *z);
839  template LIB_UTILITIES_EXPORT void Assmb(int n,
840  const Nektar::NekSingle *sign, const Nektar::NekSingle *x, const int *y,
841  Nektar::NekSingle *z);
842 
843  /************* Reduction routines *****************/
844 
845  /// \brief Subtract return sum(x)
846  template<class T> T Vsum( int n, const T *x, const int incx)
847  {
848 
849  T sum = 0;
850 
851  while( n-- )
852  {
853  sum += (*x);
854  x += incx;
855  }
856 
857  return sum;
858  }
859 
860  template LIB_UTILITIES_EXPORT Nektar::NekDouble Vsum( int n, const Nektar::NekDouble *x, const int incx);
862  const Nektar::NekSingle *x, const int incx);
863  template LIB_UTILITIES_EXPORT int Vsum( int n, const int *x, const int incx);
864 
865  /// \brief Return the index of the maximum element in x
866  template<class T> int Imax( int n, const T *x, const int incx)
867  {
868 
869  int i, indx = ( n > 0 ) ? 0 : -1;
870  T xmax = *x;
871 
872  for (i = 0; i < n; i++)
873  {
874  if (*x > xmax)
875  {
876  xmax = *x;
877  indx = i;
878  }
879  x += incx;
880  }
881 
882  return indx;
883  }
884 
885  template LIB_UTILITIES_EXPORT int Imax( int n, const Nektar::NekDouble *x, const int incx);
886  template LIB_UTILITIES_EXPORT int Imax( int n,
887  const Nektar::NekSingle *x, const int incx);
888  template LIB_UTILITIES_EXPORT int Imax( int n, const int *x, const int incx);
889 
890  /// \brief Return the maximum element in x -- called vmax to avoid
891  /// conflict with max
892  template<class T> T Vmax( int n, const T *x, const int incx)
893  {
894 
895  T xmax = *x;
896 
897  while( n-- )
898  {
899  if (*x > xmax)
900  {
901  xmax = *x;
902  }
903  x += incx;
904  }
905 
906  return xmax;
907  }
908 
909  template LIB_UTILITIES_EXPORT Nektar::NekDouble Vmax( int n, const Nektar::NekDouble *x, const int incx);
911  const Nektar::NekSingle *x, const int incx);
912  template LIB_UTILITIES_EXPORT int Vmax( int n, const int *x, const int incx);
913 
914  /// \brief Return the index of the maximum absolute element in x
915  template<class T> int Iamax( int n, const T *x, const int incx)
916  {
917 
918  int i, indx = ( n > 0 ) ? 0 : -1;
919  T xmax = *x;
920  T xm;
921 
922  for (i = 0; i < n; i++)
923  {
924  xm = (*x > 0)? *x: -*x;
925  if (xm > xmax)
926  {
927  xmax = xm;
928  indx = i;
929  }
930  x += incx;
931  }
932 
933  return indx;
934  }
935 
936  template LIB_UTILITIES_EXPORT int Iamax( int n, const Nektar::NekDouble *x, const int incx);
937  template LIB_UTILITIES_EXPORT int Iamax( int n, const Nektar::NekSingle *x,
938  const int incx);
939 
940  /// \brief Return the maximum absolute element in x
941  /// called vamax to avoid conflict with max
942  template<class T> T Vamax( int n, const T *x, const int incx)
943  {
944 
945  T xmax = *x;
946  T xm;
947 
948  while( n-- )
949  {
950  xm = (*x > 0)? *x: -*x;
951  if (xm > xmax)
952  {
953  xmax = xm;
954  }
955  x += incx;
956  }
957  return xmax;
958  }
959 
960  template LIB_UTILITIES_EXPORT Nektar::NekDouble Vamax( int n, const Nektar::NekDouble *x, const int incx);
962  const Nektar::NekSingle *x, const int incx);
963 
964 
965  /// \brief Return the index of the minimum element in x
966  template<class T> int Imin( int n, const T *x, const int incx)
967  {
968 
969  int i, indx = ( n > 0 ) ? 0 : -1;
970  T xmin = *x;
971 
972  for(i = 0;i < n;i++)
973  {
974  if( *x < xmin )
975  {
976  xmin = *x;
977  indx = i;
978  }
979  x += incx;
980  }
981 
982  return indx;
983  }
984 
985  template LIB_UTILITIES_EXPORT int Imin( int n, const Nektar::NekDouble *x, const int incx);
986  template LIB_UTILITIES_EXPORT int Imin( int n, const Nektar::NekSingle *x,
987  const int incx);
988  template LIB_UTILITIES_EXPORT int Imin( int n, const int *x, const int incx);
989 
990  /// \brief Return the minimum element in x - called vmin to avoid
991  /// conflict with min
992  template<class T> T Vmin( int n, const T *x, const int incx)
993  {
994 
995  T xmin = *x;
996 
997  while( n-- )
998  {
999  if (*x < xmin)
1000  {
1001  xmin = *x;
1002  }
1003  x += incx;
1004  }
1005 
1006  return xmin;
1007  }
1008 
1009  template LIB_UTILITIES_EXPORT Nektar::NekDouble Vmin( int n, const Nektar::NekDouble *x, const int incx);
1011  const Nektar::NekSingle *x, const int incx);
1012  template LIB_UTILITIES_EXPORT int Vmin( int n, const int *x, const int incx);
1013 
1014  /// \brief Return number of NaN elements of x
1015  template<class T> int Nnan(int n, const T *x, const int incx)
1016  {
1017 
1018  int nNan = 0;
1019 
1020  while (n--)
1021  {
1022  if (*x != *x)
1023  {
1024  nNan++;
1025  }
1026  x += incx;
1027  }
1028 
1029  return nNan;
1030  }
1031 
1032  template LIB_UTILITIES_EXPORT int Nnan(int n, const Nektar::NekDouble *x, const int incx);
1033  template LIB_UTILITIES_EXPORT int Nnan(int n, const Nektar::NekSingle *x,
1034  const int incx);
1035  template LIB_UTILITIES_EXPORT int Nnan(int n, const int *x, const int incx);
1036 
1037  /// \brief vvtvp (vector times vector times vector): z = w*x*y
1038  template<class T> T Dot( int n,
1039  const T *w,
1040  const T *x)
1041  {
1042  T sum = 0;
1043 
1044  while( n-- )
1045  {
1046  sum += (*w) * (*x);
1047  ++w;
1048  ++x;
1049  }
1050  return sum;
1051  }
1052 
1054  const Nektar::NekDouble *w,
1055  const Nektar::NekDouble *x);
1057  const Nektar::NekSingle *w,
1058  const Nektar::NekSingle *x);
1059 
1060  /// \brief vvtvp (vector times vector times vector): z = w*x*y
1061  template<class T> T Dot( int n,
1062  const T *w, const int incw,
1063  const T *x, const int incx)
1064  {
1065  T sum = 0;
1066 
1067  while( n-- )
1068  {
1069  sum += (*w) * (*x);
1070  w += incw;
1071  x += incx;
1072  }
1073  return sum;
1074  }
1075 
1077  const Nektar::NekDouble *w, const int incw,
1078  const Nektar::NekDouble *x, const int incx);
1080  const Nektar::NekSingle *w, const int incw,
1081  const Nektar::NekSingle *x, const int incx);
1082 
1083  /// \brief vvtvp (vector times vector times vector): z = w*x*y
1084  template<class T> T Dot2( int n,
1085  const T *w,
1086  const T *x,
1087  const int *y)
1088  {
1089  T sum = 0;
1090 
1091  while( n-- )
1092  {
1093  sum += (*y == 1 ? (*w) * (*x) : 0 );
1094  ++w;
1095  ++x;
1096  ++y;
1097  }
1098  return sum;
1099  }
1100 
1102  const Nektar::NekDouble *w,
1103  const Nektar::NekDouble *x,
1104  const int *y);
1106  const Nektar::NekSingle *w,
1107  const Nektar::NekSingle *x,
1108  const int *y);
1109 
1110 
1111  /// \brief vvtvp (vector times vector times vector): z = w*x*y
1112  template<class T> T Dot2( int n,
1113  const T *w, const int incw,
1114  const T *x, const int incx,
1115  const int *y, const int incy)
1116  {
1117  T sum = 0;
1118 
1119  while( n-- )
1120  {
1121  sum += (*y == 1 ? (*w) * (*x) : 0.0 );
1122  w += incw;
1123  x += incx;
1124  y += incy;
1125  }
1126  return sum;
1127  }
1128 
1130  const Nektar::NekDouble *w, const int incw,
1131  const Nektar::NekDouble *x, const int incx,
1132  const int *y, const int incy);
1134  const Nektar::NekSingle *w, const int incw,
1135  const Nektar::NekSingle *x, const int incx,
1136  const int *y, const int incy);
1137 
1138 /*
1139  // \brief copy one int vector to another
1140  void Vcopy(int n, const int *x, const int incx, int *y,
1141  const int incy)
1142  {
1143  if( incx ==1 && incy == 1)
1144  {
1145  memcpy(y,x,n*sizeof(int));
1146  }
1147  else
1148  {
1149  while( n-- )
1150  {
1151  *y = *x;
1152  x += incx;
1153  y += incy;
1154  }
1155  }
1156  }
1157 
1158  // \brief copy one unsigned int vector to another
1159  void Vcopy(int n, const unsigned int *x, const int incx, unsigned int *y,
1160  const int incy)
1161  {
1162  if( incx ==1 && incy == 1)
1163  {
1164  memcpy(y,x,n*sizeof(int));
1165  }
1166  else
1167  {
1168  while( n-- )
1169  {
1170  *y = *x;
1171  x += incx;
1172  y += incy;
1173  }
1174  }
1175  }
1176 
1177  // \brief copy one double vector to another
1178  void Vcopy(int n, const double *x, const int incx, double *y,
1179  const int incy)
1180  {
1181  if( incx ==1 && incy == 1)
1182  {
1183  memcpy(y,x,n*sizeof(double));
1184  }
1185  else
1186  {
1187  while( n-- )
1188  {
1189  *y = *x;
1190  x += incx;
1191  y += incy;
1192  }
1193  }
1194  }
1195 */
1196 
1197  // \brief copy one vector to another
1198  template<typename T>
1199  void Vcopy(int n, const T *x, const int incx,
1200  T *y, const int incy)
1201  {
1202  if( incx ==1 && incy == 1)
1203  {
1204  memcpy(y,x,n*sizeof(T));
1205  }
1206  else
1207  {
1208  while( n-- )
1209  {
1210  *y = *x;
1211  x += incx;
1212  y += incy;
1213  }
1214  }
1215  }
1216 
1217  template LIB_UTILITIES_EXPORT void Vcopy( int n, const int *x, const int incx, int *y, const int incy);
1218  template LIB_UTILITIES_EXPORT void Vcopy( int n, const unsigned int *x, const int incx, unsigned int *y, const int incy);
1219  template LIB_UTILITIES_EXPORT void Vcopy( int n, const Nektar::NekDouble *x, const int incx, Nektar::NekDouble *y, const int incy);
1220  template LIB_UTILITIES_EXPORT void Vcopy( int n,
1221  const Nektar::NekSingle *x, const int incx, Nektar::NekSingle *y,
1222  const int incy);
1223 
1224 
1225  // \brief reverse the ordering of vector to another
1226  template<class T> void Reverse( int n, const T *x, const int incx, T *y, const int incy)
1227  {
1228  int i;
1229  T store;
1230 
1231  // Perform element by element swaps in case x and y reference the same
1232  // array.
1233  int nloop = n/2;
1234 
1235  // copy value in case of n is odd number
1236  y[nloop] = x[nloop];
1237 
1238  const T* x_end = x + (n-1)*incx;
1239  T* y_end = y + (n-1)*incy;
1240  for (i = 0; i < nloop; ++i) {
1241  store = *x_end;
1242  *y_end = *x;
1243  *y = store;
1244  x += incx;
1245  y += incy;
1246  x_end -= incx;
1247  y_end -= incy;
1248  }
1249  }
1250 
1251  template LIB_UTILITIES_EXPORT void Reverse( int n, const Nektar::NekDouble *x, const int incx, Nektar::NekDouble *y, const int incy);
1252  template LIB_UTILITIES_EXPORT void Reverse( int n,
1253  const Nektar::NekSingle *x, const int incx, Nektar::NekSingle *y,
1254  const int incy);
1255 }
#define LIB_UTILITIES_EXPORT
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:15
#define NTAB
Definition: Vmath.cpp:68
#define IA2
Definition: Vmath.cpp:63
#define IR2
Definition: Vmath.cpp:67
#define NDIV
Definition: Vmath.cpp:69
#define IA1
Definition: Vmath.cpp:62
#define IM1
Definition: Vmath.cpp:58
#define IR1
Definition: Vmath.cpp:66
#define IMM1
Definition: Vmath.cpp:61
#define IQ1
Definition: Vmath.cpp:64
#define RNMX
Definition: Vmath.cpp:71
#define IQ2
Definition: Vmath.cpp:65
#define AM
Definition: Vmath.cpp:60
#define IM2
Definition: Vmath.cpp:59
double NekDouble
Definition: Vmath.cpp:40
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:475
void Ssub(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha - x.
Definition: Vmath.cpp:405
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:691
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:192
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:565
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
Definition: Vmath.cpp:493
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:1084
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:461
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:992
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:513
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:846
void Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
Definition: Vmath.cpp:772
T Dot(int n, const T *w, const T *x)
vvtvp (vector times vector times vector): z = w*x*y
Definition: Vmath.cpp:1038
void Gathr(int n, const T *sign, const T *x, const int *y, T *z)
Gather vector z[i] = sign[i]*x[y[i]].
Definition: Vmath.cpp:756
void Assmb(int n, const T *x, const int *y, T *z)
Assemble z[y[i]] += x[i]; z should be zero'd first.
Definition: Vmath.cpp:813
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:602
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:322
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:541
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:658
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:225
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:291
int Imax(int n, const T *x, const int incx)
Return the index of the maximum element in x.
Definition: Vmath.cpp:866
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:257
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
Definition: Vmath.cpp:966
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:436
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
T ran2(long *idum)
Generates a number from ~Normal(0,1)
Definition: Vmath.cpp:74
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:723
void FillWhiteNoise(int n, const T eps, T *x, const int incx, int outseed)
Fills a vector with white noise.
Definition: Vmath.cpp:142
int Nnan(int n, const T *x, const int incx)
Return number of NaN elements of x.
Definition: Vmath.cpp:1015
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:942
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:341
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1226
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:892
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:625
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1199
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:372
int Iamax(int n, const T *x, const int incx)
Return the index of the maximum absolute element in x.
Definition: Vmath.cpp:915
scalarT< T > log(scalarT< T > in)
Definition: scalar.hpp:278
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:267