Nektar++
VmathArray.hpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File VmathArray.hpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Wrappers around Vmath routines using Array<OneD,T> as arugments
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #ifndef NEKTAR_LIB_LIBUTILITIES_BASSICUTILS_VECTORMATHARRAY_HPP
36 #define NEKTAR_LIB_LIBUTILITIES_BASSICUTILS_VECTORMATHARRAY_HPP
37 
38 #include <LibUtilities/BasicUtils/SharedArray.hpp>
39 #include <LibUtilities/BasicUtils/Vmath.hpp>
40 #include <LibUtilities/BasicUtils/VmathSIMD.hpp>
41 
42  namespace Vmath
43  {
44  using namespace Nektar;
45 
46  /***************** Math routines ***************/
47  /// \brief Fill a vector with a constant value
48  template<class T> void Fill( int n, const T alpha, Array<OneD, T> &x, const int incx )
49  {
50 
51  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Out of bounds");
52 
53  Fill(n,alpha,&x[0],incx);
54  }
55 
56  template<class T> void FillWhiteNoise(
57  int n, const T eps, Array<OneD, T> &x,
58  const int incx, int outseed = 9999)
59  {
60  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Out of bounds");
61 
62  FillWhiteNoise(n,eps,&x[0],incx,outseed);
63  }
64 
65  /// \brief Multiply vector z = x*y
66  template<class T> void Vmul( int n, const Array<OneD, const T> &x, const int incx, const Array<OneD, const T> &y, const int incy, Array<OneD,T> &z, const int incz)
67  {
68  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Array out of bounds");
69  ASSERTL1(n*incy <= y.size()+y.GetOffset(),"Array out of bounds");
70  ASSERTL1(n*incz <= z.size()+z.GetOffset(),"Array out of bounds");
71 
72 #ifdef NEKTAR_ENABLE_SIMD_VMATH
73  boost::ignore_unused(incx, incy, incz);
74  ASSERTL1(incx == 1, "Simd vmath requires inc = 1");
75  ASSERTL1(incy == 1, "Simd vmath requires inc = 1");
76  ASSERTL1(incz == 1, "Simd vmath requires inc = 1");
77  SIMD::Vmul(n,&x[0],&y[0],&z[0]);
78 #else
79  Vmul(n,&x[0],incx,&y[0],incy,&z[0],incz);
80 #endif
81 
82  }
83 
84  template<class T> void Vmul( int n, const Array<TwoD,NekDouble>::const_reference &x, const int incx, const Array<OneD,const T> &y, const int incy, Array<OneD,T> &z, const int incz)
85  {
86  ASSERTL1(n*incx <= x.size(),"Array out of bounds");
87  ASSERTL1(n*incy <= y.size()+y.GetOffset(),"Array out of bounds");
88  ASSERTL1(n*incz <= z.size()+z.GetOffset(),"Array out of bounds");
89 
90  Vmul(n,x.origin(),incx,&y[0],incy,&z[0],incz);
91  }
92 
93  /// \brief Scalar multiply y = alpha*y
94 
95  template<class T> void Smul( int n, const T alpha, const Array<OneD,const T> &x, const int incx, Array<OneD,T> &y, const int incy)
96  {
97  ASSERTL1(static_cast<unsigned int>(n*incx) <= x.size()+x.GetOffset(),"Array out of bounds");
98  ASSERTL1(static_cast<unsigned int>(n*incy) <= y.size()+y.GetOffset(),"Array out of bounds");
99 
100  Smul(n,alpha, &x[0],incx,&y[0],incy);
101  }
102 
103  /// \brief Multiply vector z = x/y
104  template<class T> void Vdiv( int n, const Array<OneD,const T> &x, const int incx, const Array<OneD,const T> &y, const int incy, Array<OneD,T> &z, const int incz)
105  {
106  ASSERTL1(static_cast<unsigned int>(n*incx) <= x.size()+x.GetOffset(),"Array out of bounds");
107  ASSERTL1(static_cast<unsigned int>(n*incy) <= y.size()+y.GetOffset(),"Array out of bounds");
108  ASSERTL1(static_cast<unsigned int>(n*incz) <= z.size()+z.GetOffset(),"Array out of bounds");
109 
110  Vdiv(n,&x[0],incx,&y[0],incy,&z[0],incz);
111 
112  }
113 
114  /// \brief Scalar multiply y = alpha/y
115  template<class T> void Sdiv( int n, const T alpha, const Array<OneD,const T> &x, const int incx, Array<OneD,T> &y, const int incy)
116  {
117  ASSERTL1(static_cast<unsigned int>(n*incx) <= x.size()+x.GetOffset(),"Array out of bounds");
118  ASSERTL1(static_cast<unsigned int>(n*incy) <= y.size()+y.GetOffset(),"Array out of bounds");
119 
120  Sdiv(n,alpha,&x[0],incx,&y[0],incy);
121  }
122 
123  /// \brief Add vector z = x+y
124  template<class T> void Vadd( int n, const Array<OneD,const T> &x, const int incx, const Array<OneD,const T> &y, const int incy, Array<OneD,T> &z, const int incz)
125  {
126  ASSERTL1(static_cast<unsigned int>(n*incx) <= x.size()+x.GetOffset(),"Array out of bounds");
127  ASSERTL1(static_cast<unsigned int>(n*incy) <= y.size()+y.GetOffset(),"Array out of bounds");
128  ASSERTL1(static_cast<unsigned int>(n*incz) <= z.size()+z.GetOffset(),"Array out of bounds");
129 
130 #ifdef NEKTAR_ENABLE_SIMD_VMATH
131  boost::ignore_unused(incx, incy, incz);
132  ASSERTL1(incx == 1, "Simd vmath requires inc = 1");
133  ASSERTL1(incy == 1, "Simd vmath requires inc = 1");
134  ASSERTL1(incz == 1, "Simd vmath requires inc = 1");
135  SIMD::Vadd(n,&x[0],&y[0],&z[0]);
136 #else
137  Vadd(n,&x[0],incx,&y[0],incy,&z[0],incz);
138 #endif
139  }
140 
141  /// \brief Add vector y = alpha + x
142  template<class T> void Sadd( int n, const T alpha, const Array<OneD,const T> &x,const int incx, Array<OneD,T> &y, const int incy)
143  {
144 
145  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Array out of bounds");
146  ASSERTL1(n*incy <= y.size()+y.GetOffset(),"Array out of bounds");
147 
148  Sadd(n,alpha,&x[0],incx,&y[0],incy);
149  }
150 
151  /// \brief Subtract vector z = x-y
152  template<class T> void Vsub( int n, const Array<OneD,const T> &x, const int incx, const Array<OneD,const T> &y, const int incy, Array<OneD,T> &z, const int incz)
153  {
154  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Array out of bounds");
155  ASSERTL1(n*incy <= y.size()+y.GetOffset(),"Array out of bounds");
156  ASSERTL1(n*incz <= z.size()+z.GetOffset(),"Array out of bounds");
157 
158  Vsub(n,&x[0],incx,&y[0],incy,&z[0],incz);
159 
160  }
161 
162  /// \brief Add vector y = alpha - x
163  template<class T> void Ssub( int n, const T alpha,
164  const Array<OneD,const T> &x,const int incx, Array<OneD,T> &y,
165  const int incy)
166  {
167 
168  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Array out of bounds");
169  ASSERTL1(n*incy <= y.size()+y.GetOffset(),"Array out of bounds");
170 
171  Ssub(n,alpha,&x[0],incx,&y[0],incy);
172  }
173 
174  /// \brief Zero vector
175  template<class T> void Zero(int n, Array<OneD,T> &x, const int incx)
176  {
177  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Array out of bounds");
178 
179  Zero(n,&x[0],incx);
180 
181  }
182 
183  /// \brief Negate x = -x
184  template<class T> void Neg( int n, Array<OneD,T> &x, const int incx)
185  {
186  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Array out of bounds");
187 
188  Neg(n,&x[0],incx);
189 
190  }
191 
192  template<class T> void Vlog(int n, const Array<OneD,const T> &x, const int incx, Array<OneD,T> &y, const int incy)
193  {
194  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Array out of bounds");
195  ASSERTL1(n*incy <= y.size()+y.GetOffset(),"Array out of bounds");
196 
197  Vlog(n, &x[0], incx, &y[0], incy);
198  }
199 
200 
201  template<class T> void Vexp(int n, const Array<OneD,const T> &x, const int incx, Array<OneD,T> &y, const int incy)
202  {
203  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Array out of bounds");
204  ASSERTL1(n*incy <= y.size()+y.GetOffset(),"Array out of bounds");
205 
206  Vexp(n, &x[0], incx, &y[0], incy);
207  }
208 
209  template<class T> void Vpow(int n, const Array<OneD,const T> &x, const int incx, const T f, Array<OneD,T> &y, const int incy)
210  {
211  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Array out of bounds");
212  ASSERTL1(n*incy <= y.size()+y.GetOffset(),"Array out of bounds");
213 
214  Vpow(n, &x[0], incx, f, &y[0], incy);
215  }
216 
217  /// \brief sqrt y = sqrt(x)
218  template<class T> void Vsqrt(int n, const Array<OneD,const T> &x, const int incx, Array<OneD,T> &y, const int incy)
219  {
220  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Array out of bounds");
221  ASSERTL1(n*incy <= y.size()+y.GetOffset(),"Array out of bounds");
222 
223  Vsqrt(n,&x[0],incx,&y[0],incy);
224  }
225 
226  /// \brief vabs: y = |x|
227  template<class T> void Vabs(int n, const Array<OneD,const T> &x, const int incx, Array<OneD,T> &y, const int incy)
228  {
229  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Array out of bounds");
230  ASSERTL1(n*incy <= y.size()+y.GetOffset(),"Array out of bounds");
231 
232  Vabs(n,&x[0],incx,&y[0],incy);
233  }
234 
235  /********** Triad routines ***********************/
236 
237  /// \brief vvtvp (vector times vector plus vector): z = w*x + y
238  template<class T> void Vvtvp(int n, const Array<OneD, const T> &w, const int incw, const Array<OneD,const T> &x, const int incx, const Array<OneD, const T> &y, const int incy, Array<OneD,T> &z, const int incz)
239  {
240  ASSERTL1(n*incw <= w.size()+w.GetOffset(),"Array out of bounds");
241  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Array out of bounds");
242  ASSERTL1(n*incy <= y.size()+y.GetOffset(),"Array out of bounds");
243  ASSERTL1(n*incz <= z.size()+z.GetOffset(),"Array out of bounds");
244 
245 #ifdef NEKTAR_ENABLE_SIMD_VMATH
246  boost::ignore_unused(incw, incx, incy, incz);
247  ASSERTL1(incw == 1, "Simd vmath requires inc = 1");
248  ASSERTL1(incx == 1, "Simd vmath requires inc = 1");
249  ASSERTL1(incy == 1, "Simd vmath requires inc = 1");
250  ASSERTL1(incz == 1, "Simd vmath requires inc = 1");
251  SIMD::Vvtvp(n,&w[0],&x[0],&y[0],&z[0]);
252 #else
253  Vvtvp(n,&w[0],incw,&x[0],incx,&y[0],incy,&z[0],incz);
254 #endif
255  }
256 
257  template<class T> void Vvtvp(int n, const Array<TwoD,NekDouble>::const_reference &w, const int incw, const Array<OneD, const T> &x, const int incx, const Array<OneD,const T> &y, const int incy, Array<OneD,T> &z, const int incz)
258  {
259  ASSERTL1(n*incw <= w.size(),"Array out of bounds");
260  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Array out of bounds");
261  ASSERTL1(n*incy <= y.size()+y.GetOffset(),"Array out of bounds");
262  ASSERTL1(n*incz <= z.size()+z.GetOffset(),"Array out of bounds");
263 
264  Vvtvp(n,w.origin(),incw,&x[0],incx,&y[0],incy,&z[0],incz);
265  }
266 
267  /// \brief svtvp (scalar times vector plus vector): z = alpha*x + y
268  template<class T> void Svtvp(int n, const T alpha, const Array<OneD,const T> &x, const int incx, const Array<OneD, const T> &y, const int incy, Array<OneD,T> &z, const int incz)
269  {
270  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Array out of bounds");
271  ASSERTL1(n*incy <= y.size()+y.GetOffset(),"Array out of bounds");
272  ASSERTL1(n*incz <= z.size()+z.GetOffset(),"Array out of bounds");
273 
274  Svtvp(n,alpha,&x[0],incx,&y[0],incy,&z[0],incz);
275 
276  }
277 
278  /// \brief svtvp (scalar times vector plus vector): z = alpha*x + y
279  template<class T> void Svtvm(int n, const T alpha, const Array<OneD,const T> &x, const int incx, const Array<OneD, const T> &y, const int incy, Array<OneD,T> &z, const int incz)
280  {
281  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Array out of bounds");
282  ASSERTL1(n*incy <= y.size()+y.GetOffset(),"Array out of bounds");
283  ASSERTL1(n*incz <= z.size()+z.GetOffset(),"Array out of bounds");
284 
285  Svtvm(n,alpha,&x[0],incx,&y[0],incy,&z[0],incz);
286 
287  }
288 
289  /// \brief vvtvm (vector times vector minus vector): z = w*x - y
290  template<class T> void Vvtvm(int n, const Array<OneD,const T> &w, const int incw, const Array<OneD,const T> &x, const int incx, const Array<OneD,const T> &y, const int incy, Array<OneD,T> &z, const int incz)
291  {
292  ASSERTL1(n*incw <= w.size()+w.GetOffset(),"Array out of bounds");
293  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Array out of bounds");
294  ASSERTL1(n*incy <= y.size()+y.GetOffset(),"Array out of bounds");
295  ASSERTL1(n*incz <= z.size()+z.GetOffset(),"Array out of bounds");
296 
297 
298 
299 #ifdef NEKTAR_ENABLE_SIMD_VMATH
300  boost::ignore_unused(incw, incx, incy, incz);
301  ASSERTL1(incw == 1, "Simd vmath requires inc = 1");
302  ASSERTL1(incx == 1, "Simd vmath requires inc = 1");
303  ASSERTL1(incy == 1, "Simd vmath requires inc = 1");
304  ASSERTL1(incz == 1, "Simd vmath requires inc = 1");
305  SIMD::Vvtvm(n,&w[0],&x[0],&y[0],&z[0]);
306 #else
307  Vvtvm(n,&w[0],incw,&x[0],incx,&y[0],incy,&z[0],incz);
308 #endif
309 
310  }
311 
312  /// \brief vvtvvtp (vector times vector plus vector times vector): z = v*w + y*z
313  template<class T> void Vvtvvtp (
314  int n,
315  const Array<OneD,const T> &v, int incv,
316  const Array<OneD,const T> &w, int incw,
317  const Array<OneD,const T> &x, int incx,
318  const Array<OneD,const T> &y, int incy,
319  Array<OneD, T> &z, int incz)
320  {
321  ASSERTL1(n*incv <= v.size()+v.GetOffset(),"Array out of bounds");
322  ASSERTL1(n*incw <= w.size()+w.GetOffset(),"Array out of bounds");
323  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Array out of bounds");
324  ASSERTL1(n*incy <= y.size()+y.GetOffset(),"Array out of bounds");
325  ASSERTL1(n*incz <= z.size()+z.GetOffset(),"Array out of bounds");
326 
327 #ifdef NEKTAR_ENABLE_SIMD_VMATH
328  boost::ignore_unused(incv, incw, incx, incy, incz);
329  ASSERTL1(incw == 1, "Simd vmath requires inc = 1");
330  ASSERTL1(incx == 1, "Simd vmath requires inc = 1");
331  ASSERTL1(incy == 1, "Simd vmath requires inc = 1");
332  ASSERTL1(incz == 1, "Simd vmath requires inc = 1");
333  SIMD::Vvtvvtp(n,&v[0],&w[0],&x[0],&y[0],&z[0]);
334 #else
335  Vvtvvtp(n,&v[0],incv,&w[0],incw,&x[0],incx,&y[0],incy,&z[0],incz);
336 #endif
337  }
338 
339  /// \brief svtsvtp (scalar times vector plus scalar times vector): z = alpha*x + beta*y
340  template<class T> void Svtsvtp(int n, const T alpha, const Array<OneD,const T> &x, const int incx, const T beta, const Array<OneD,const T> &y, const int incy, Array<OneD,T> &z, const int incz)
341  {
342  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Array out of bounds");
343  ASSERTL1(n*incy <= y.size()+y.GetOffset(),"Array out of bounds");
344  ASSERTL1(n*incz <= z.size()+z.GetOffset(),"Array out of bounds");
345 
346  Svtsvtp(n,alpha,&x[0],incx,beta,&y[0],incy,&z[0],incz);
347  }
348 
349 
350 
351 
352  /************ Misc routine from Veclib (and extras) ************/
353 
354  /// \brief Gather vector z[i] = x[y[i]]
355  template<class T, class I, typename = typename std::enable_if
356  <
357  std::is_floating_point<T>::value &&
358  std::is_integral<I>::value
359  >::type
360  >
361  void Gathr(I n, const Array<OneD, const T> &x, const Array<OneD, I> &y,
362  Array<OneD,T> &z)
363  {
364  ASSERTL1(n <= y.size()+y.GetOffset(),"Array out of bounds");
365  ASSERTL1(n <= z.size()+z.GetOffset(),"Array out of bounds");
366 
367 #ifdef NEKTAR_ENABLE_SIMD_VMATH
368  SIMD::Gathr(n,&x[0],&y[0],&z[0]);
369 #else
370  Gathr(n,&x[0],&y[0],&z[0]);
371 #endif
372  }
373 
374  /// \brief Scatter vector z[y[i]] = x[i]
375  template<class T> void Scatr(int n, const Array<OneD,const T> &x, const Array<OneD,const int> &y, Array<OneD,T> &z)
376  {
377  ASSERTL1(n <= x.size()+x.GetOffset(),"Array out of bounds");
378  ASSERTL1(n <= y.size()+y.GetOffset(),"Array out of bounds");
379 
380  Scatr(n,&x[0],&y[0],&z[0]);
381  }
382 
383 
384  /// \brief Assemble z[y[i]] += x[i]; z should be zero'd first
385  template<class T> void Assmb(int n, const Array<OneD,T> &x, const Array<OneD,int> &y, Array<OneD,T> &z)
386  {
387  ASSERTL1(n <= x.size()+x.GetOffset(),"Array out of bounds");
388  ASSERTL1(n <= y.size()+y.GetOffset(),"Array out of bounds");
389 
390  Assmb(n,&x[0],&y[0],&z[0]);
391  }
392 
393 
394  /************* Reduction routines *****************/
395 
396  /// \brief Subtract return sum(x)
397  template<class T> T Vsum( int n, const Array<OneD, const T> &x, const int incx)
398  {
399  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Array out of bounds");
400 
401  return Vsum(n,&x[0],incx);
402  }
403 
404 
405  /// \brief Return the index of the maximum element in x
406  template<class T> int Imax( int n, const Array<OneD, const T> &x, const int incx)
407  {
408  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Array out of bounds");
409 
410  return Imax(n,&x[0],incx);
411  }
412 
413  /// \brief Return the maximum element in x -- called vmax to avoid
414  /// conflict with max
415  template<class T> T Vmax( int n, const Array<OneD, const T> &x, const int incx)
416  {
417  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Array out of bounds");
418 
419  return Vmax(n,&x[0],incx);
420  }
421 
422  /// \brief Return the index of the maximum absolute element in x
423  template<class T> int Iamax( int n, const Array<OneD, const T> &x, const int incx)
424  {
425  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Array out of bounds");
426 
427  return Iamax(n,&x[0],incx);
428 
429  }
430 
431  /// \brief Return the maximum absolute element in x
432  /// called vamax to avoid conflict with max
433  template<class T> T Vamax( int n, const Array<OneD, const T> &x, const int incx)
434  {
435  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Array out of bounds");
436 
437  return Vamax(n,&x[0],incx);
438  }
439 
440 
441  /// \brief Return the index of the minimum element in x
442  template<class T> int Imin( int n, const Array<OneD, const T> &x, const int incx)
443  {
444  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Array out of bounds");
445 
446  return Imin(n,&x[0],incx);
447  }
448 
449  /// \brief Return the minimum element in x - called vmin to avoid
450  /// conflict with min
451  template<class T> T Vmin( int n, const Array<OneD, const T> &x, const int incx)
452  {
453  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Array out of bounds");
454 
455  return Vmin(n,&x[0],incx);
456  }
457 
458  /// \brief Return number of NaN elements of x
459  template<class T> int Nnan(int n, const Array<OneD, const T> &x, const int incx)
460  {
461  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
462 
463  return Nnan(n, &x[0], incx);
464  }
465 
466  /// \brief
467  template<class T> T Dot(int n,
468  const Array<OneD, const T> &w,
469  const Array<OneD, const T> &x)
470  {
471  ASSERTL1(n <= w.size()+w.GetOffset(),"Array out of bounds");
472  ASSERTL1(n <= x.size()+x.GetOffset(),"Array out of bounds");
473 
474  return Dot(n,&w[0],&x[0]);
475  }
476 
477  /// \brief
478  template<class T> T Dot(int n,
479  const Array<OneD, const T> &w, const int incw,
480  const Array<OneD, const T> &x, const int incx)
481  {
482  ASSERTL1(n*incw <= w.size()+w.GetOffset(),"Array out of bounds");
483  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Array out of bounds");
484 
485  return Dot(n,&w[0],incw,&x[0],incx);
486  }
487 
488  /// \brief
489  template<class T> T Dot2(int n,
490  const Array<OneD, const T> &w,
491  const Array<OneD, const T> &x,
492  const Array<OneD, const int> &y)
493  {
494  ASSERTL1(n <= w.size()+w.GetOffset(),"Array out of bounds");
495  ASSERTL1(n <= x.size()+x.GetOffset(),"Array out of bounds");
496  ASSERTL1(n <= y.size()+y.GetOffset(),"Array out of bounds");
497 
498  return Dot2(n,&w[0],&x[0],&y[0]);
499  }
500 
501  /// \brief
502  template<class T> T Ddot(int n,
503  const Array<OneD, const T> &w, const int incw,
504  const Array<OneD, const T> &x, const int incx,
505  const Array<OneD, const int> &y, const int incy)
506  {
507  ASSERTL1(n*incw <= w.size()+w.GetOffset(),"Array out of bounds");
508  ASSERTL1(n*incx <= x.size()+x.GetOffset(),"Array out of bounds");
509  ASSERTL1(n*incy <= y.size()+y.GetOffset(),"Array out of bounds");
510 
511  return Dot2(n,&w[0],incw,&x[0],incx,&y[0],incy);
512  }
513 
514  /********** Memory routines ***********************/
515 
516  template<class T> void Vcopy(int n, const Array<OneD, const T> &x, int incx, Array<OneD,T> &y, int const incy)
517  {
518  ASSERTL1(static_cast<unsigned int>(std::abs(n*incx)) <= x.size()+x.GetOffset(),"Array out of bounds");
519  ASSERTL1(static_cast<unsigned int>(std::abs(n*incy)) <= y.size()+y.GetOffset(),"Array out of bounds");
520 
521  Vcopy(n,&x[0],incx,&y[0],incy);
522  }
523 
524  template<class T> void Reverse(int n, const Array<OneD, const T> &x, int incx, Array<OneD,T> &y, int const incy)
525  {
526  ASSERTL1(static_cast<unsigned int>(std::abs(n*incx)) <= x.size()+x.GetOffset(),"Array out of bounds");
527  ASSERTL1(static_cast<unsigned int>(std::abs(n*incy)) <= y.size()+y.GetOffset(),"Array out of bounds");
528 
529  Reverse(n,&x[0],incx,&y[0],incy);
530  }
531 
532  }
533 #endif //VECTORMATHARRAY_HPP
534 
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:250
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
Vmath::SIMD::Vvtvp
void Vvtvp(const size_t n, const T *w, const T *x, const T *y, T *z)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: VmathSIMD.hpp:264
Vmath::SIMD::Vadd
void Vadd(const size_t n, const T *x, const T *y, T *z)
Multiply vector z = x + y.
Definition: VmathSIMD.hpp:52
Vmath::SIMD::Vvtvm
void Vvtvm(const size_t n, const T *w, const T *x, const T *y, T *z)
vvtvm (vector times vector plus vector): z = w*x - y
Definition: VmathSIMD.hpp:312
Vmath::SIMD::Gathr
void Gathr(const I n, const T *x, const I *y, T *z)
Gather vector z[i] = x[y[i]].
Definition: VmathSIMD.hpp:419
Vmath::SIMD::Vvtvvtp
void Vvtvvtp(const size_t n, const T *v, const T *w, const T *x, const T *y, T *z)
vvtvvtp (vector times vector plus vector times vector):
Definition: VmathSIMD.hpp:361
Vmath::SIMD::Vmul
void Vmul(const size_t n, const T *x, const T *y, T *z)
Multiply vector z = x * y.
Definition: VmathSIMD.hpp:158
Vmath
Definition: Vmath.cpp:40
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:475
Vmath::Ssub
void Ssub(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha - x.
Definition: Vmath.cpp:405
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:691
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:192
Vmath::Vlog
void Vlog(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:104
Vmath::Vexp
void Vexp(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:116
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:565
Vmath::Vabs
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
Definition: Vmath.cpp:493
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:1084
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:461
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:992
Vmath::Ddot
T Ddot(int n, const Array< OneD, const T > &w, const int incw, const Array< OneD, const T > &x, const int incx, const Array< OneD, const int > &y, const int incy)
Definition: VmathArray.hpp:502
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:513
Vmath::Vsum
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:846
Vmath::Scatr
void Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
Definition: Vmath.cpp:772
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:1038
Vmath::Gathr
void Gathr(int n, const T *sign, const T *x, const int *y, T *z)
Gather vector z[i] = sign[i]*x[y[i]].
Definition: Vmath.cpp:756
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:813
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:602
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:322
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:541
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:225
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:291
Vmath::Imax
int Imax(int n, const T *x, const int incx)
Return the index of the maximum element in x.
Definition: Vmath.cpp:866
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:257
Vmath::Imin
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
Definition: Vmath.cpp:966
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:436
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
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:142
Vmath::Nnan
int Nnan(int n, const T *x, const int incx)
Return number of NaN elements of x.
Definition: Vmath.cpp:1015
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:942
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:341
Vmath::Reverse
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1226
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:892
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:625
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1199
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:372
Vmath::Vpow
void Vpow(int n, const T *x, const int incx, const T f, T *y, const int incy)
Definition: Vmath.hpp:127
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:915
tinysimd::abs
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:272
Generated on Wed Nov 24 2021 10:33:29 for Nektar++ by doxygen 1.9.1