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