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