Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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,
139  const int incx, int outseed)
140  {
141  // Protect the static vars here and in ran2
142  boost::mutex::scoped_lock l(mutex);
143 
144  // Define static variables for generating random numbers
145  static int iset = 0;
146  static T gset;
147  static long seed = 0;
148 
149  // Bypass seed if outseed was specified
150  if( outseed != 9999)
151  {
152  seed = long(outseed);
153  }
154 
155  while( n-- )
156  {
157  T fac, rsq, v1, v2;
158 
159  if (iset == 0)
160  {
161  do
162  {
163  v1 = 2.0 * ran2<T> (&seed) - 1.0;
164  v2 = 2.0 * ran2<T> (&seed) - 1.0;
165  rsq = v1*v1 + v2*v2;
166  } while (rsq >= 1.0 || rsq == 0.0);
167  fac = sqrt(-2.0 * log (rsq) / rsq);
168  gset = v1 * fac;
169  iset = 1;
170  *x = eps * v2 * fac;
171  }
172  else
173  {
174  iset = 0;
175  *x = eps * gset;
176  }
177  x += incx;
178  }
179  }
180  template LIB_UTILITIES_EXPORT void FillWhiteNoise( int n, const Nektar::NekDouble eps, Nektar::NekDouble *x, const int incx, int outseed);
181 
182  /// \brief Multiply vector z = x*y
183  template<class T> void Vmul( int n, const T *x, const int incx, const T *y,
184  const int incy, T*z, const int incz)
185  {
186  ++n;
187  if (incx == 1 && incy == 1 && incz == 1)
188  {
189  while( --n )
190  {
191  *z = (*x) * (*y);
192  ++x;
193  ++y;
194  ++z;
195  }
196  }
197  else
198  {
199  while( --n )
200  {
201  *z = (*x) * (*y);
202  x += incx;
203  y += incy;
204  z += incz;
205  }
206  }
207  }
208  template LIB_UTILITIES_EXPORT void Vmul( int n, const Nektar::NekDouble *x, const int incx, const Nektar::NekDouble *y,
209  const int incy, Nektar::NekDouble*z, const int incz);
210 
211  /// \brief Scalar multiply y = alpha*y
212 
213  template<class T> void Smul( int n, const T alpha, const T *x, const int incx,
214  T *y, const int incy)
215  {
216  ++n;
217  if (incx == 1 && incy == 1)
218  {
219  while( --n )
220  {
221  *y = alpha * (*x);
222  ++x;
223  ++y;
224  }
225  }
226  else
227  {
228  while( --n )
229  {
230  *y = alpha * (*x);
231  x += incx;
232  y += incy;
233  }
234  }
235  }
236 
237  template LIB_UTILITIES_EXPORT void Smul( int n, const Nektar::NekDouble alpha, const Nektar::NekDouble *x, const int incx,
238  Nektar::NekDouble *y, const int incy);
239 
240  /// \brief Multiply vector z = x/y
241  template<class T> void Vdiv( int n, const T *x, const int incx, const T *y,
242  const int incy, T*z, const int incz)
243  {
244  ++n;
245  if (incx == 1 && incy == 1)
246  {
247  while( --n )
248  {
249  *z = (*x) / (*y);
250  ++x;
251  ++y;
252  ++z;
253  }
254  }
255  else
256  {
257  while( --n )
258  {
259  *z = (*x) / (*y);
260  x += incx;
261  y += incy;
262  z += incz;
263  }
264  }
265  }
266 
267  template LIB_UTILITIES_EXPORT void Vdiv( int n, const Nektar::NekDouble *x, const int incx, const Nektar::NekDouble *y,
268  const int incy, Nektar::NekDouble*z, const int incz);
269 
270  /// \brief Scalar multiply y = alpha/y
271  template<class T> void Sdiv( int n, const T alpha, const T *x,
272  const int incx, T *y, const int incy)
273  {
274  ++n;
275  if (incx == 1 && incy == 1)
276  {
277  while( --n )
278  {
279  *y = alpha / (*x);
280  ++x;
281  ++y;
282  }
283  }
284  else
285  {
286  while( --n )
287  {
288  *y = alpha / (*x);
289  x += incx;
290  y += incy;
291  }
292  }
293  }
294 
295  template LIB_UTILITIES_EXPORT void Sdiv( int n, const Nektar::NekDouble alpha, const Nektar::NekDouble *x,
296  const int incx, Nektar::NekDouble *y, const int incy);
297 
298  /// \brief Add vector z = x+y
299  template<class T> void Vadd( int n, const T *x, const int incx, const T *y,
300  const int incy, T *z, const int incz)
301  {
302  while( n-- )
303  {
304  *z = (*x) + (*y);
305  x += incx;
306  y += incy;
307  z += incz;
308  }
309  }
310 
311  template LIB_UTILITIES_EXPORT void Vadd( int n, const Nektar::NekDouble *x, const int incx, const Nektar::NekDouble *y,
312  const int incy, Nektar::NekDouble *z, const int incz);
313 
314  /// \brief Add vector y = alpha + x
315  template<class T> void Sadd( int n, const T alpha, const T *x,
316  const int incx, T *y, const int incy)
317  {
318  ++n;
319  if (incx == 1 && incy == 1)
320  {
321  while( --n )
322  {
323  *y = alpha + (*x);
324  ++x;
325  ++y;
326  }
327  }
328  else
329  {
330  while( --n )
331  {
332  *y = alpha + (*x);
333  x += incx;
334  y += incy;
335  }
336  }
337  }
338 
339  template LIB_UTILITIES_EXPORT void Sadd( int n, const Nektar::NekDouble alpha, const Nektar::NekDouble *x,
340  const int incx, Nektar::NekDouble *y, const int incy);
341 
342  /// \brief Subtract vector z = x-y
343  template<class T> void Vsub( int n, const T *x, const int incx, const T *y,
344  const int incy, T *z, const int incz)
345  {
346  ++n;
347  if (incx == 1 && incy == 1 && incz == 1)
348  {
349  while( --n )
350  {
351  *z = (*x) - (*y);
352  ++x;
353  ++y;
354  ++z;
355  }
356  }
357  else
358  {
359  while( --n )
360  {
361  *z = (*x) - (*y);
362  x += incx;
363  y += incy;
364  z += incz;
365  }
366  }
367  }
368 
369  template LIB_UTILITIES_EXPORT void Vsub( int n, const Nektar::NekDouble *x, const int incx, const Nektar::NekDouble *y,
370  const int incy, Nektar::NekDouble *z, const int incz);
371 
372  /// \brief Zero vector
373  template<class T> void Zero(int n, T *x, const int incx)
374  {
375  if(incx == 1)
376  {
377  std::memset(x,'\0', n*sizeof(T));
378  }
379  else
380  {
381  T zero = 0;
382  ++n;
383  while(--n)
384  {
385  *x = zero;
386  x+=incx;
387  }
388  }
389  }
390 
391  template LIB_UTILITIES_EXPORT void Zero(int n, Nektar::NekDouble *x, const int incx);
392  template LIB_UTILITIES_EXPORT void Zero(int n, int *x, const int incx);
393  template LIB_UTILITIES_EXPORT void Zero(int n, long *x, const int incx);
394 
395  /// \brief Negate x = -x
396  template<class T> void Neg( int n, T *x, const int incx)
397  {
398  while( n-- )
399  {
400  *x = -(*x);
401  x += incx;
402  }
403  }
404 
405  template LIB_UTILITIES_EXPORT void Neg( int n, Nektar::NekDouble *x, const int incx);
406 
407  /// \brief sqrt y = sqrt(x)
408  template<class T> void Vsqrt(int n, const T *x, const int incx,
409  T *y, const int incy)
410  {
411  while (n--)
412  {
413  *y = sqrt( *x );
414  x += incx;
415  y += incy;
416  }
417  }
418 
419  template LIB_UTILITIES_EXPORT void Vsqrt(int n, const Nektar::NekDouble *x, const int incx,
420  Nektar::NekDouble *y, const int incy);
421 
422 
423  /// \brief vabs: y = |x|
424  template<class T> void Vabs(int n, const T *x, const int incx,
425  T *y, const int incy)
426  {
427  while( n-- )
428  {
429  *y = ( *x >0)? *x:-(*x);
430  x += incx;
431  y += incy;
432  }
433  }
434 
435  template LIB_UTILITIES_EXPORT void Vabs(int n, const Nektar::NekDouble *x, const int incx,
436  Nektar::NekDouble *y, const int incy);
437 
438 
439  /********** Triad routines ***********************/
440 
441  /// \brief vvtvp (vector times vector plus vector): z = w*x + y
442  template<class T> void Vvtvp(int n,
443  const T *w, const int incw,
444  const T *x, const int incx,
445  const T *y, const int incy,
446  T *z, const int incz)
447  {
448  while( n-- )
449  {
450  *z = (*w) * (*x) + (*y);
451  w += incw;
452  x += incx;
453  y += incy;
454  z += incz;
455  }
456  }
457 
458  template LIB_UTILITIES_EXPORT void Vvtvp(int n,
459  const Nektar::NekDouble *w, const int incw,
460  const Nektar::NekDouble *x, const int incx,
461  const Nektar::NekDouble *y, const int incy,
462  Nektar::NekDouble *z, const int incz);
463 
464  /// \brief vvtvm (vector times vector plus vector): z = w*x - y
465  template<class T> void Vvtvm(int n, const T *w, const int incw, const T *x,
466  const int incx, const T *y, const int incy,
467  T *z, const int incz)
468  {
469  while( n-- )
470  {
471  *z = (*w) * (*x) - (*y);
472  w += incw;
473  x += incx;
474  y += incy;
475  z += incz;
476  }
477  }
478 
479  template LIB_UTILITIES_EXPORT void Vvtvm(int n, const Nektar::NekDouble *w, const int incw, const Nektar::NekDouble *x,
480  const int incx, const Nektar::NekDouble *y, const int incy,
481  Nektar::NekDouble *z, const int incz);
482 
483 
484  /// \brief svtvp (scalar times vector plus vector): z = alpha*x + y
485  template<class T> LIB_UTILITIES_EXPORT void Svtvp(int n, const T alpha, const T *x,
486  const int incx, const T *y, const int incy,
487  T *z, const int incz)
488  {
489  ++n;
490  if (incx == 1 && incy == 1 && incz == 1)
491  {
492  while( --n )
493  {
494  *z = alpha * (*x) + (*y);
495  ++x;
496  ++y;
497  ++z;
498  }
499  }
500  else
501  {
502  while( --n )
503  {
504  *z = alpha * (*x) + (*y);
505  x += incx;
506  y += incy;
507  z += incz;
508  }
509  }
510  }
511 
512  template LIB_UTILITIES_EXPORT void Svtvp(int n, const Nektar::NekDouble alpha, const Nektar::NekDouble *x,
513  const int incx, const Nektar::NekDouble *y, const int incy,
514  Nektar::NekDouble *z, const int incz);
515 
516 
517  /// \brief svtvp (scalar times vector plus vector): z = alpha*x - y
518  template<class T> void Svtvm(int n, const T alpha, const T *x,
519  const int incx, const T *y, const int incy,
520  T *z, const int incz)
521  {
522  while( n-- )
523  {
524  *z = alpha * (*x) - (*y);
525  x += incx;
526  y += incy;
527  z += incz;
528  }
529  }
530 
531  template LIB_UTILITIES_EXPORT void Svtvm(int n, const Nektar::NekDouble alpha, const Nektar::NekDouble *x,
532  const int incx, const Nektar::NekDouble *y, const int incy,
533  Nektar::NekDouble *z, const int incz);
534 
535  /// \brief vvtvvtp (vector times vector plus vector times vector):
536  // z = v*w + x*y
537  template<class T> void Vvtvvtp (int n,
538  const T* v, int incv,
539  const T* w, int incw,
540  const T* x, int incx,
541  const T* y, int incy,
542  T* z, int incz)
543  {
544  while( n-- )
545  {
546  *z = (*v) * (*w) + (*x) * (*y);
547  v += incv;
548  w += incw;
549  x += incx;
550  y += incy;
551  z += incz;
552  }
553  }
554  template LIB_UTILITIES_EXPORT void Vvtvvtp (int n,
555  const Nektar::NekDouble* v, int incv,
556  const Nektar::NekDouble* w, int incw,
557  const Nektar::NekDouble* x, int incx,
558  const Nektar::NekDouble* y, int incy,
559  Nektar::NekDouble* z, int incz);
560 
561 
562  /// \brief vvtvvtm (vector times vector minus vector times vector):
563  // z = v*w - x*y
564  template<class T> void Vvtvvtm (int n,
565  const T* v, int incv,
566  const T* w, int incw,
567  const T* x, int incx,
568  const T* y, int incy,
569  T* z, int incz)
570  {
571  while( n-- )
572  {
573  *z = (*v) * (*w) - (*x) * (*y);
574  v += incv;
575  w += incw;
576  x += incx;
577  y += incy;
578  z += incz;
579  }
580  }
581 
582  template LIB_UTILITIES_EXPORT void Vvtvvtm (int n,
583  const Nektar::NekDouble* v, int incv,
584  const Nektar::NekDouble* w, int incw,
585  const Nektar::NekDouble* x, int incx,
586  const Nektar::NekDouble* y, int incy,
587  Nektar::NekDouble* z, int incz);
588 
589  /// \brief vvtvvtp (scalar times vector plus scalar times vector):
590  // z = alpha*x + beta*y
591  template<class T> void Svtsvtp (int n,
592  const T alpha,
593  const T* x, int incx,
594  const T beta,
595  const T* y, int incy,
596  T* z, int incz)
597  {
598  while( n-- )
599  {
600  *z = alpha * (*x) + beta * (*y);
601  x += incx;
602  y += incy;
603  z += incz;
604  }
605  }
606 
607  template LIB_UTILITIES_EXPORT void Svtsvtp (int n,
608  const Nektar::NekDouble alpha,
609  const Nektar::NekDouble* x, int incx,
610  const Nektar::NekDouble beta,
611  const Nektar::NekDouble* y, int incy,
612  Nektar::NekDouble* z, int incz);
613 
614 
615  /// \brief Vstvpp (scalar times vector plus vector plus vector):
616  // z = v*w + x*y
617  template<class T> void Vstvpp(int n,
618  const T alpha,
619  const T* v, int incv,
620  const T* w, int incw,
621  const T* x, int incx,
622  T* z, int incz)
623  {
624  while( n-- )
625  {
626  *z = alpha * (*v) + (*w) + (*x);
627  v += incv;
628  w += incw;
629  x += incx;
630  z += incz;
631  }
632  }
633 
634  template LIB_UTILITIES_EXPORT void Vstvpp(int n,
635  const Nektar::NekDouble alpha,
636  const Nektar::NekDouble* v, int incv,
637  const Nektar::NekDouble* w, int incw,
638  const Nektar::NekDouble* x, int incx,
639  Nektar::NekDouble* z, int incz);
640 
641  /************ Misc routine from Veclib (and extras) ************/
642 
643  /// \brief Gather vector z[i] = x[y[i]]
644  template<class T> void Gathr(int n, const T *x, const int *y,
645  T *z)
646  {
647  while (n--)
648  {
649  *z++ = *(x + *y++);
650  }
651  return;
652  }
653 
654  template LIB_UTILITIES_EXPORT void Gathr(int n, const Nektar::NekDouble *x, const int *y,
655  Nektar::NekDouble *z);
656 
657 
658  /// \brief Gather vector z[i] = sign[i]*x[y[i]]
659  template<class T> void Gathr(int n, const T *sign, const T *x, const int *y,
660  T *z)
661  {
662  while (n--)
663  {
664  *z++ = *(sign++) * (*(x + *y++));
665  }
666  return;
667  }
668 
669  template LIB_UTILITIES_EXPORT void Gathr(int n, const Nektar::NekDouble *sign, const Nektar::NekDouble *x, const int *y,
670  Nektar::NekDouble *z);
671 
672  /// \brief Scatter vector z[y[i]] = x[i]
673  template<class T> void Scatr(int n, const T *x, const int *y,
674  T *z)
675  {
676  while (n--)
677  {
678  *(z + *(y++)) = *(x++);
679  }
680  }
681 
682  template LIB_UTILITIES_EXPORT void Scatr(int n, const Nektar::NekDouble *x, const int *y,
683  Nektar::NekDouble *z);
684 
685  /// \brief Scatter vector z[y[i]] = sign[i]*x[i]
686  template<class T> void Scatr(int n, const T *sign, const T *x, const int *y,
687  T *z)
688  {
689  while (n--)
690  {
691  if(*sign)
692  {
693  *(z + *(y++)) = *(sign++) * (*(x++));
694  }
695  else
696  {
697  x++;
698  y++;
699  sign++;
700  }
701  }
702  }
703 
704  template LIB_UTILITIES_EXPORT void Scatr(int n, const Nektar::NekDouble *sign, const Nektar::NekDouble *x, const int *y,
705  Nektar::NekDouble *z);
706 
707 
708  /// \brief Assemble z[y[i]] += x[i]; z should be zero'd first
709  template<class T> void Assmb(int n, const T *x, const int *y,
710  T *z)
711  {
712  while (n--)
713  {
714  *(z + *(y++)) += *(x++);
715  }
716  }
717 
718  template LIB_UTILITIES_EXPORT void Assmb(int n, const Nektar::NekDouble *x, const int *y,
719  Nektar::NekDouble *z);
720 
721  /// \brief Assemble z[y[i]] += sign[i]*x[i]; z should be zero'd first
722  template<class T> void Assmb(int n, const T *sign, const T *x, const int *y,
723  T *z)
724  {
725  while (n--)
726  {
727  *(z + *(y++)) += *(sign++) * (*(x++));
728  }
729  }
730 
731  template LIB_UTILITIES_EXPORT void Assmb(int n, const Nektar::NekDouble *sign, const Nektar::NekDouble *x, const int *y,
732  Nektar::NekDouble *z);
733 
734  /************* Reduction routines *****************/
735 
736  /// \brief Subtract return sum(x)
737  template<class T> T Vsum( int n, const T *x, const int incx)
738  {
739 
740  T sum = 0;
741 
742  while( n-- )
743  {
744  sum += (*x);
745  x += incx;
746  }
747 
748  return sum;
749  }
750 
751  template LIB_UTILITIES_EXPORT Nektar::NekDouble Vsum( int n, const Nektar::NekDouble *x, const int incx);
752  template LIB_UTILITIES_EXPORT int Vsum( int n, const int *x, const int incx);
753 
754  /// \brief Return the index of the maximum element in x
755  template<class T> int Imax( int n, const T *x, const int incx)
756  {
757 
758  int i, indx = ( n > 0 ) ? 0 : -1;
759  T xmax = *x;
760 
761  for (i = 0; i < n; i++)
762  {
763  if (*x > xmax)
764  {
765  xmax = *x;
766  indx = i;
767  }
768  x += incx;
769  }
770 
771  return indx;
772  }
773 
774  template LIB_UTILITIES_EXPORT int Imax( int n, const Nektar::NekDouble *x, const int incx);
775  template LIB_UTILITIES_EXPORT int Imax( int n, const int *x, const int incx);
776 
777  /// \brief Return the maximum element in x -- called vmax to avoid
778  /// conflict with max
779  template<class T> T Vmax( int n, const T *x, const int incx)
780  {
781 
782  T xmax = *x;
783 
784  while( n-- )
785  {
786  if (*x > xmax)
787  {
788  xmax = *x;
789  }
790  x += incx;
791  }
792 
793  return xmax;
794  }
795 
796  template LIB_UTILITIES_EXPORT Nektar::NekDouble Vmax( int n, const Nektar::NekDouble *x, const int incx);
797  template LIB_UTILITIES_EXPORT int Vmax( int n, const int *x, const int incx);
798 
799  /// \brief Return the index of the maximum absolute element in x
800  template<class T> int Iamax( int n, const T *x, const int incx)
801  {
802 
803  int i, indx = ( n > 0 ) ? 0 : -1;
804  T xmax = *x;
805  T xm;
806 
807  for (i = 0; i < n; i++)
808  {
809  xm = (*x > 0)? *x: -*x;
810  if (xm > xmax)
811  {
812  xmax = xm;
813  indx = i;
814  }
815  x += incx;
816  }
817 
818  return indx;
819  }
820 
821  template LIB_UTILITIES_EXPORT int Iamax( int n, const Nektar::NekDouble *x, const int incx);
822 
823  /// \brief Return the maximum absolute element in x
824  /// called vamax to avoid conflict with max
825  template<class T> T Vamax( int n, const T *x, const int incx)
826  {
827 
828  T xmax = *x;
829  T xm;
830 
831  while( n-- )
832  {
833  xm = (*x > 0)? *x: -*x;
834  if (xm > xmax)
835  {
836  xmax = xm;
837  }
838  x += incx;
839  }
840  return xmax;
841  }
842 
843  template LIB_UTILITIES_EXPORT Nektar::NekDouble Vamax( int n, const Nektar::NekDouble *x, const int incx);
844 
845 
846  /// \brief Return the index of the minimum element in x
847  template<class T> int Imin( int n, const T *x, const int incx)
848  {
849 
850  int i, indx = ( n > 0 ) ? 0 : -1;
851  T xmin = *x;
852 
853  for(i = 0;i < n;i++)
854  {
855  if( *x < xmin )
856  {
857  xmin = *x;
858  indx = i;
859  }
860  x += incx;
861  }
862 
863  return indx;
864  }
865 
866  template LIB_UTILITIES_EXPORT int Imin( int n, const Nektar::NekDouble *x, const int incx);
867  template LIB_UTILITIES_EXPORT int Imin( int n, const int *x, const int incx);
868 
869  /// \brief Return the minimum element in x - called vmin to avoid
870  /// conflict with min
871  template<class T> T Vmin( int n, const T *x, const int incx)
872  {
873 
874  T xmin = *x;
875 
876  while( n-- )
877  {
878  if (*x < xmin)
879  {
880  xmin = *x;
881  }
882  x += incx;
883  }
884 
885  return xmin;
886  }
887 
888  template LIB_UTILITIES_EXPORT Nektar::NekDouble Vmin( int n, const Nektar::NekDouble *x, const int incx);
889  template LIB_UTILITIES_EXPORT int Vmin( int n, const int *x, const int incx);
890 
891  /// \brief Return number of NaN elements of x
892  template<class T> int Nnan(int n, const T *x, const int incx)
893  {
894 
895  int nNan = 0;
896 
897  while (n--)
898  {
899  if (*x != *x)
900  {
901  nNan++;
902  }
903  x += incx;
904  }
905 
906  return nNan;
907  }
908 
909  template LIB_UTILITIES_EXPORT int Nnan(int n, const Nektar::NekDouble *x, const int incx);
910  template LIB_UTILITIES_EXPORT int Nnan(int n, const float *x, const int incx);
911  template LIB_UTILITIES_EXPORT int Nnan(int n, const int *x, const int incx);
912 
913  /// \brief vvtvp (vector times vector times vector): z = w*x*y
914  template<class T> T Dot( int n,
915  const T *w,
916  const T *x)
917  {
918  T sum = 0;
919 
920  while( n-- )
921  {
922  sum += (*w) * (*x);
923  ++w;
924  ++x;
925  }
926  return sum;
927  }
928 
929  template LIB_UTILITIES_EXPORT Nektar::NekDouble Dot( int n,
930  const Nektar::NekDouble *w,
931  const Nektar::NekDouble *x);
932 
933  /// \brief vvtvp (vector times vector times vector): z = w*x*y
934  template<class T> T Dot( int n,
935  const T *w, const int incw,
936  const T *x, const int incx)
937  {
938  T sum = 0;
939 
940  while( n-- )
941  {
942  sum += (*w) * (*x);
943  w += incw;
944  x += incx;
945  }
946  return sum;
947  }
948 
949  template LIB_UTILITIES_EXPORT Nektar::NekDouble Dot( int n,
950  const Nektar::NekDouble *w, const int incw,
951  const Nektar::NekDouble *x, const int incx);
952 
953  /// \brief vvtvp (vector times vector times vector): z = w*x*y
954  template<class T> T Dot2( int n,
955  const T *w,
956  const T *x,
957  const int *y)
958  {
959  T sum = 0;
960 
961  while( n-- )
962  {
963  sum += (*y == 1 ? (*w) * (*x) : 0 );
964  ++w;
965  ++x;
966  ++y;
967  }
968  return sum;
969  }
970 
971  template LIB_UTILITIES_EXPORT Nektar::NekDouble Dot2( int n,
972  const Nektar::NekDouble *w,
973  const Nektar::NekDouble *x,
974  const int *y);
975 
976 
977  /// \brief vvtvp (vector times vector times vector): z = w*x*y
978  template<class T> T Dot2( int n,
979  const T *w, const int incw,
980  const T *x, const int incx,
981  const int *y, const int incy)
982  {
983  T sum = 0;
984 
985  while( n-- )
986  {
987  sum += (*y == 1 ? (*w) * (*x) : 0.0 );
988  w += incw;
989  x += incx;
990  y += incy;
991  }
992  return sum;
993  }
994 
995  template LIB_UTILITIES_EXPORT Nektar::NekDouble Dot2( int n,
996  const Nektar::NekDouble *w, const int incw,
997  const Nektar::NekDouble *x, const int incx,
998  const int *y, const int incy);
999 
1000 /*
1001  // \brief copy one int vector to another
1002  void Vcopy(int n, const int *x, const int incx, int *y,
1003  const int incy)
1004  {
1005  if( incx ==1 && incy == 1)
1006  {
1007  memcpy(y,x,n*sizeof(int));
1008  }
1009  else
1010  {
1011  while( n-- )
1012  {
1013  *y = *x;
1014  x += incx;
1015  y += incy;
1016  }
1017  }
1018  }
1019 
1020  // \brief copy one unsigned int vector to another
1021  void Vcopy(int n, const unsigned int *x, const int incx, unsigned int *y,
1022  const int incy)
1023  {
1024  if( incx ==1 && incy == 1)
1025  {
1026  memcpy(y,x,n*sizeof(int));
1027  }
1028  else
1029  {
1030  while( n-- )
1031  {
1032  *y = *x;
1033  x += incx;
1034  y += incy;
1035  }
1036  }
1037  }
1038 
1039  // \brief copy one double vector to another
1040  void Vcopy(int n, const double *x, const int incx, double *y,
1041  const int incy)
1042  {
1043  if( incx ==1 && incy == 1)
1044  {
1045  memcpy(y,x,n*sizeof(double));
1046  }
1047  else
1048  {
1049  while( n-- )
1050  {
1051  *y = *x;
1052  x += incx;
1053  y += incy;
1054  }
1055  }
1056  }
1057 */
1058 
1059  // \brief copy one vector to another
1060  template<typename T>
1061  void Vcopy(int n, const T *x, const int incx,
1062  T *y, const int incy)
1063  {
1064  if( incx ==1 && incy == 1)
1065  {
1066  memcpy(y,x,n*sizeof(T));
1067  }
1068  else
1069  {
1070  while( n-- )
1071  {
1072  *y = *x;
1073  x += incx;
1074  y += incy;
1075  }
1076  }
1077  }
1078 
1079  template LIB_UTILITIES_EXPORT void Vcopy( int n, const int *x, const int incx, int *y, const int incy);
1080  template LIB_UTILITIES_EXPORT void Vcopy( int n, const unsigned int *x, const int incx, unsigned int *y, const int incy);
1081  template LIB_UTILITIES_EXPORT void Vcopy( int n, const Nektar::NekDouble *x, const int incx, Nektar::NekDouble *y, const int incy);
1082 
1083 
1084  // \brief reverse the ordering of vector to another
1085  template<class T> void Reverse( int n, const T *x, const int incx, T *y, const int incy)
1086  {
1087  int i;
1088  T store;
1089  int nloop = n/2;
1090 
1091  // copy value in case of n is odd number
1092  y[nloop] = x[nloop];
1093 
1094  for(i = 0; i < nloop; ++i)
1095  {
1096  store = x[n-1-i];
1097  y[n-1-i] = x[i];
1098  y[i] = store;
1099  }
1100  }
1101 
1102  template LIB_UTILITIES_EXPORT void Reverse( int n, const Nektar::NekDouble *x, const int incx, Nektar::NekDouble *y, const int incy);
1103 }
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:644
sign
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:27
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:408
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:779
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:871
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:847
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:485
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:442
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:271
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:241
Vmath::Vabs
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
Definition: Vmath.cpp:424
Vmath::Imax
int Imax(int n, const T *x, const int incx)
Return the index of the maximum element in x.
Definition: Vmath.cpp:755
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:1085
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:213
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:892
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:518
Vmath::Scatr
void Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
Definition: Vmath.cpp:673
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:396
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:709
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:564
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:825
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:800
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:315
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:914
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:343
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:617
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:537
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:465
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:954
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:591
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:737
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
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:299
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:183
Generated on Sat Nov 25 2017 16:18:19 for Nektar++ by   doxygen 1.8.8