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 // 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  *(z + *(y++)) = *(sign++) * (*(x++));
678  }
679  }
680 
681  template LIB_UTILITIES_EXPORT void Scatr(int n, const Nektar::NekDouble *sign, const Nektar::NekDouble *x, const int *y,
682  Nektar::NekDouble *z);
683 
684 
685  /// \brief Assemble z[y[i]] += x[i]; z should be zero'd first
686  template<class T> void Assmb(int n, const T *x, const int *y,
687  T *z)
688  {
689  while (n--)
690  {
691  *(z + *(y++)) += *(x++);
692  }
693  }
694 
695  template LIB_UTILITIES_EXPORT void Assmb(int n, const Nektar::NekDouble *x, const int *y,
696  Nektar::NekDouble *z);
697 
698  /// \brief Assemble z[y[i]] += sign[i]*x[i]; z should be zero'd first
699  template<class T> void Assmb(int n, const T *sign, const T *x, const int *y,
700  T *z)
701  {
702  while (n--)
703  {
704  *(z + *(y++)) += *(sign++) * (*(x++));
705  }
706  }
707 
708  template LIB_UTILITIES_EXPORT void Assmb(int n, const Nektar::NekDouble *sign, const Nektar::NekDouble *x, const int *y,
709  Nektar::NekDouble *z);
710 
711  /************* Reduction routines *****************/
712 
713  /// \brief Subtract return sum(x)
714  template<class T> T Vsum( int n, const T *x, const int incx)
715  {
716 
717  T sum = 0;
718 
719  while( n-- )
720  {
721  sum += (*x);
722  x += incx;
723  }
724 
725  return sum;
726  }
727 
728  template LIB_UTILITIES_EXPORT Nektar::NekDouble Vsum( int n, const Nektar::NekDouble *x, const int incx);
729  template LIB_UTILITIES_EXPORT int Vsum( int n, const int *x, const int incx);
730 
731  /// \brief Return the index of the maximum element in x
732  template<class T> int Imax( int n, const T *x, const int incx)
733  {
734 
735  int i, indx = ( n > 0 ) ? 0 : -1;
736  T xmax = *x;
737 
738  for (i = 0; i < n; i++)
739  {
740  if (*x > xmax)
741  {
742  xmax = *x;
743  indx = i;
744  }
745  x += incx;
746  }
747 
748  return indx;
749  }
750 
751  template LIB_UTILITIES_EXPORT int Imax( int n, const Nektar::NekDouble *x, const int incx);
752  template LIB_UTILITIES_EXPORT int Imax( int n, const int *x, const int incx);
753 
754  /// \brief Return the maximum element in x -- called vmax to avoid
755  /// conflict with max
756  template<class T> T Vmax( int n, const T *x, const int incx)
757  {
758 
759  T xmax = *x;
760 
761  while( n-- )
762  {
763  if (*x > xmax)
764  {
765  xmax = *x;
766  }
767  x += incx;
768  }
769 
770  return xmax;
771  }
772 
773  template LIB_UTILITIES_EXPORT Nektar::NekDouble Vmax( int n, const Nektar::NekDouble *x, const int incx);
774  template LIB_UTILITIES_EXPORT int Vmax( int n, const int *x, const int incx);
775 
776  /// \brief Return the index of the maximum absolute element in x
777  template<class T> int Iamax( int n, const T *x, const int incx)
778  {
779 
780  int i, indx = ( n > 0 ) ? 0 : -1;
781  T xmax = *x;
782  T xm;
783 
784  for (i = 0; i < n; i++)
785  {
786  xm = (*x > 0)? *x: -*x;
787  if (xm > xmax)
788  {
789  xmax = xm;
790  indx = i;
791  }
792  x += incx;
793  }
794 
795  return indx;
796  }
797 
798  template LIB_UTILITIES_EXPORT int Iamax( int n, const Nektar::NekDouble *x, const int incx);
799 
800  /// \brief Return the maximum absolute element in x
801  /// called vamax to avoid conflict with max
802  template<class T> T Vamax( int n, const T *x, const int incx)
803  {
804 
805  T xmax = *x;
806  T xm;
807 
808  while( n-- )
809  {
810  xm = (*x > 0)? *x: -*x;
811  if (xm > xmax)
812  {
813  xmax = xm;
814  }
815  x += incx;
816  }
817  return xmax;
818  }
819 
820  template LIB_UTILITIES_EXPORT Nektar::NekDouble Vamax( int n, const Nektar::NekDouble *x, const int incx);
821 
822 
823  /// \brief Return the index of the minimum element in x
824  template<class T> int Imin( int n, const T *x, const int incx)
825  {
826 
827  int i, indx = ( n > 0 ) ? 0 : -1;
828  T xmin = *x;
829 
830  for(i = 0;i < n;i++)
831  {
832  if( *x < xmin )
833  {
834  xmin = *x;
835  indx = i;
836  }
837  x += incx;
838  }
839 
840  return indx;
841  }
842 
843  template LIB_UTILITIES_EXPORT int Imin( int n, const Nektar::NekDouble *x, const int incx);
844  template LIB_UTILITIES_EXPORT int Imin( int n, const int *x, const int incx);
845 
846  /// \brief Return the minimum element in x - called vmin to avoid
847  /// conflict with min
848  template<class T> T Vmin( int n, const T *x, const int incx)
849  {
850 
851  T xmin = *x;
852 
853  while( n-- )
854  {
855  if (*x < xmin)
856  {
857  xmin = *x;
858  }
859  x += incx;
860  }
861 
862  return xmin;
863  }
864 
865  template LIB_UTILITIES_EXPORT Nektar::NekDouble Vmin( int n, const Nektar::NekDouble *x, const int incx);
866  template LIB_UTILITIES_EXPORT int Vmin( int n, const int *x, const int incx);
867 
868  /// \brief Return number of NaN elements of x
869  template<class T> int Nnan(int n, const T *x, const int incx)
870  {
871 
872  int nNan = 0;
873 
874  while (n--)
875  {
876  if (*x != *x)
877  {
878  nNan++;
879  }
880  x += incx;
881  }
882 
883  return nNan;
884  }
885 
886  template LIB_UTILITIES_EXPORT int Nnan(int n, const Nektar::NekDouble *x, const int incx);
887  template LIB_UTILITIES_EXPORT int Nnan(int n, const float *x, const int incx);
888  template LIB_UTILITIES_EXPORT int Nnan(int n, const int *x, const int incx);
889 
890  /// \brief vvtvp (vector times vector times vector): z = w*x*y
891  template<class T> T Dot( int n,
892  const T *w,
893  const T *x)
894  {
895  T sum = 0;
896 
897  while( n-- )
898  {
899  sum += (*w) * (*x);
900  ++w;
901  ++x;
902  }
903  return sum;
904  }
905 
906  template LIB_UTILITIES_EXPORT Nektar::NekDouble Dot( int n,
907  const Nektar::NekDouble *w,
908  const Nektar::NekDouble *x);
909 
910  /// \brief vvtvp (vector times vector times vector): z = w*x*y
911  template<class T> T Dot( int n,
912  const T *w, const int incw,
913  const T *x, const int incx)
914  {
915  T sum = 0;
916 
917  while( n-- )
918  {
919  sum += (*w) * (*x);
920  w += incw;
921  x += incx;
922  }
923  return sum;
924  }
925 
926  template LIB_UTILITIES_EXPORT Nektar::NekDouble Dot( int n,
927  const Nektar::NekDouble *w, const int incw,
928  const Nektar::NekDouble *x, const int incx);
929 
930  /// \brief vvtvp (vector times vector times vector): z = w*x*y
931  template<class T> T Dot2( int n,
932  const T *w,
933  const T *x,
934  const int *y)
935  {
936  T sum = 0;
937 
938  while( n-- )
939  {
940  sum += (*y == 1 ? (*w) * (*x) : 0 );
941  ++w;
942  ++x;
943  ++y;
944  }
945  return sum;
946  }
947 
948  template LIB_UTILITIES_EXPORT Nektar::NekDouble Dot2( int n,
949  const Nektar::NekDouble *w,
950  const Nektar::NekDouble *x,
951  const int *y);
952 
953 
954  /// \brief vvtvp (vector times vector times vector): z = w*x*y
955  template<class T> T Dot2( int n,
956  const T *w, const int incw,
957  const T *x, const int incx,
958  const int *y, const int incy)
959  {
960  T sum = 0;
961 
962  while( n-- )
963  {
964  sum += (*y == 1 ? (*w) * (*x) : 0.0 );
965  w += incw;
966  x += incx;
967  y += incy;
968  }
969  return sum;
970  }
971 
972  template LIB_UTILITIES_EXPORT Nektar::NekDouble Dot2( int n,
973  const Nektar::NekDouble *w, const int incw,
974  const Nektar::NekDouble *x, const int incx,
975  const int *y, const int incy);
976 
977 /*
978  // \brief copy one int vector to another
979  void Vcopy(int n, const int *x, const int incx, int *y,
980  const int incy)
981  {
982  if( incx ==1 && incy == 1)
983  {
984  memcpy(y,x,n*sizeof(int));
985  }
986  else
987  {
988  while( n-- )
989  {
990  *y = *x;
991  x += incx;
992  y += incy;
993  }
994  }
995  }
996 
997  // \brief copy one unsigned int vector to another
998  void Vcopy(int n, const unsigned int *x, const int incx, unsigned int *y,
999  const int incy)
1000  {
1001  if( incx ==1 && incy == 1)
1002  {
1003  memcpy(y,x,n*sizeof(int));
1004  }
1005  else
1006  {
1007  while( n-- )
1008  {
1009  *y = *x;
1010  x += incx;
1011  y += incy;
1012  }
1013  }
1014  }
1015 
1016  // \brief copy one double vector to another
1017  void Vcopy(int n, const double *x, const int incx, double *y,
1018  const int incy)
1019  {
1020  if( incx ==1 && incy == 1)
1021  {
1022  memcpy(y,x,n*sizeof(double));
1023  }
1024  else
1025  {
1026  while( n-- )
1027  {
1028  *y = *x;
1029  x += incx;
1030  y += incy;
1031  }
1032  }
1033  }
1034 */
1035 
1036  // \brief copy one vector to another
1037  template<typename T>
1038  void Vcopy(int n, const T *x, const int incx,
1039  T *y, const int incy)
1040  {
1041  if( incx ==1 && incy == 1)
1042  {
1043  memcpy(y,x,n*sizeof(T));
1044  }
1045  else
1046  {
1047  while( n-- )
1048  {
1049  *y = *x;
1050  x += incx;
1051  y += incy;
1052  }
1053  }
1054  }
1055 
1056  template LIB_UTILITIES_EXPORT void Vcopy( int n, const int *x, const int incx, int *y, const int incy);
1057  template LIB_UTILITIES_EXPORT void Vcopy( int n, const unsigned int *x, const int incx, unsigned int *y, const int incy);
1058  template LIB_UTILITIES_EXPORT void Vcopy( int n, const Nektar::NekDouble *x, const int incx, Nektar::NekDouble *y, const int incy);
1059 
1060 
1061  // \brief reverse the ordering of vector to another
1062  template<class T> void Reverse( int n, const T *x, const int incx, T *y, const int incy)
1063  {
1064  int i;
1065  T store;
1066  int nloop = n/2;
1067 
1068  // copy value in case of n is odd number
1069  y[nloop] = x[nloop];
1070 
1071  for(i = 0; i < nloop; ++i)
1072  {
1073  store = x[n-1-i];
1074  y[n-1-i] = x[i];
1075  y[i] = store;
1076  }
1077  }
1078 
1079  template LIB_UTILITIES_EXPORT void Reverse( int n, const Nektar::NekDouble *x, const int incx, Nektar::NekDouble *y, const int incy);
1080 }
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:756
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:848
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:824
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:732
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:1062
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:869
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:43
Vmath::Assmb
void Assmb(int n, const T *x, const int *y, T *z)
Assemble z[y[i]] += x[i]; z should be zero'd first.
Definition: Vmath.cpp:686
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:802
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:777
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:891
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:931
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:714
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:1038
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 Tue Dec 8 2015 23:45:31 for Nektar++ by   doxygen 1.8.9.1