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>
49 void Fill(int n, const T alpha, Array<OneD, T> &x, const int incx)
50 {
51 
52  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Out of bounds");
53 
54  Fill(n, alpha, &x[0], incx);
55 }
56 
57 template <class T>
58 void FillWhiteNoise(int n, const T eps, Array<OneD, T> &x, const int incx,
59  int outseed = 9999)
60 {
61  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Out of bounds");
62 
63  FillWhiteNoise(n, eps, &x[0], incx, outseed);
64 }
65 
66 /// \brief Multiply vector z = x*y
67 template <class T>
68 void Vmul(int n, const Array<OneD, const T> &x, const int incx,
69  const Array<OneD, const T> &y, const int incy, Array<OneD, T> &z,
70  const int incz)
71 {
72  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
73  ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
74  ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
75 
76 #ifdef NEKTAR_ENABLE_SIMD_VMATH
77  boost::ignore_unused(incx, incy, incz);
78  ASSERTL1(incx == 1, "Simd vmath requires inc = 1");
79  ASSERTL1(incy == 1, "Simd vmath requires inc = 1");
80  ASSERTL1(incz == 1, "Simd vmath requires inc = 1");
81  SIMD::Vmul(n, &x[0], &y[0], &z[0]);
82 #else
83  Vmul(n, &x[0], incx, &y[0], incy, &z[0], incz);
84 #endif
85 }
86 
87 template <class T>
88 void Vmul(int n, const Array<TwoD, NekDouble>::const_reference &x,
89  const int incx, const Array<OneD, const T> &y, const int incy,
90  Array<OneD, T> &z, const int incz)
91 {
92  ASSERTL1(n * incx <= x.size(), "Array out of bounds");
93  ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
94  ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
95 
96  Vmul(n, x.origin(), incx, &y[0], incy, &z[0], incz);
97 }
98 
99 /// \brief Scalar multiply y = alpha*y
100 
101 template <class T>
102 void Smul(int n, const T alpha, const Array<OneD, const T> &x, const int incx,
103  Array<OneD, T> &y, const int incy)
104 {
105  ASSERTL1(static_cast<unsigned int>(n * incx) <= x.size() + x.GetOffset(),
106  "Array out of bounds");
107  ASSERTL1(static_cast<unsigned int>(n * incy) <= y.size() + y.GetOffset(),
108  "Array out of bounds");
109 
110  Smul(n, alpha, &x[0], incx, &y[0], incy);
111 }
112 
113 /// \brief Multiply vector z = x/y
114 template <class T>
115 void Vdiv(int n, const Array<OneD, const T> &x, const int incx,
116  const Array<OneD, const T> &y, const int incy, Array<OneD, T> &z,
117  const int incz)
118 {
119  ASSERTL1(static_cast<unsigned int>(n * incx) <= x.size() + x.GetOffset(),
120  "Array out of bounds");
121  ASSERTL1(static_cast<unsigned int>(n * incy) <= y.size() + y.GetOffset(),
122  "Array out of bounds");
123  ASSERTL1(static_cast<unsigned int>(n * incz) <= z.size() + z.GetOffset(),
124  "Array out of bounds");
125 
126  Vdiv(n, &x[0], incx, &y[0], incy, &z[0], incz);
127 }
128 
129 /// \brief Scalar multiply y = alpha/y
130 template <class T>
131 void Sdiv(int n, const T alpha, const Array<OneD, const T> &x, const int incx,
132  Array<OneD, T> &y, const int incy)
133 {
134  ASSERTL1(static_cast<unsigned int>(n * incx) <= x.size() + x.GetOffset(),
135  "Array out of bounds");
136  ASSERTL1(static_cast<unsigned int>(n * incy) <= y.size() + y.GetOffset(),
137  "Array out of bounds");
138 
139  Sdiv(n, alpha, &x[0], incx, &y[0], incy);
140 }
141 
142 /// \brief Add vector z = x+y
143 template <class T>
144 void Vadd(int n, const Array<OneD, const T> &x, const int incx,
145  const Array<OneD, const T> &y, const int incy, Array<OneD, T> &z,
146  const int incz)
147 {
148  ASSERTL1(static_cast<unsigned int>(n * incx) <= x.size() + x.GetOffset(),
149  "Array out of bounds");
150  ASSERTL1(static_cast<unsigned int>(n * incy) <= y.size() + y.GetOffset(),
151  "Array out of bounds");
152  ASSERTL1(static_cast<unsigned int>(n * incz) <= z.size() + z.GetOffset(),
153  "Array out of bounds");
154 
155 #ifdef NEKTAR_ENABLE_SIMD_VMATH
156  boost::ignore_unused(incx, incy, incz);
157  ASSERTL1(incx == 1, "Simd vmath requires inc = 1");
158  ASSERTL1(incy == 1, "Simd vmath requires inc = 1");
159  ASSERTL1(incz == 1, "Simd vmath requires inc = 1");
160  SIMD::Vadd(n, &x[0], &y[0], &z[0]);
161 #else
162  Vadd(n, &x[0], incx, &y[0], incy, &z[0], incz);
163 #endif
164 }
165 
166 /// \brief Add vector y = alpha + x
167 template <class T>
168 void Sadd(int n, const T alpha, const Array<OneD, const T> &x, const int incx,
169  Array<OneD, T> &y, const int incy)
170 {
171 
172  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
173  ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
174 
175  Sadd(n, alpha, &x[0], incx, &y[0], incy);
176 }
177 
178 /// \brief Subtract vector z = x-y
179 template <class T>
180 void Vsub(int n, const Array<OneD, const T> &x, const int incx,
181  const Array<OneD, const T> &y, const int incy, Array<OneD, T> &z,
182  const int incz)
183 {
184  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
185  ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
186  ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
187 
188  Vsub(n, &x[0], incx, &y[0], incy, &z[0], incz);
189 }
190 
191 /// \brief Add vector y = alpha - x
192 template <class T>
193 void Ssub(int n, const T alpha, const Array<OneD, const T> &x, const int incx,
194  Array<OneD, T> &y, const int incy)
195 {
196 
197  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
198  ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
199 
200  Ssub(n, alpha, &x[0], incx, &y[0], incy);
201 }
202 
203 /// \brief Zero vector
204 template <class T> void Zero(int n, Array<OneD, T> &x, const int incx)
205 {
206  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
207 
208  Zero(n, &x[0], incx);
209 }
210 
211 /// \brief Negate x = -x
212 template <class T> void Neg(int n, Array<OneD, T> &x, const int incx)
213 {
214  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
215 
216  Neg(n, &x[0], incx);
217 }
218 
219 template <class T>
220 void Vlog(int n, const Array<OneD, const T> &x, const int incx,
221  Array<OneD, T> &y, const int incy)
222 {
223  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
224  ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
225 
226  Vlog(n, &x[0], incx, &y[0], incy);
227 }
228 
229 template <class T>
230 void Vexp(int n, const Array<OneD, const T> &x, const int incx,
231  Array<OneD, T> &y, const int incy)
232 {
233  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
234  ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
235 
236  Vexp(n, &x[0], incx, &y[0], incy);
237 }
238 
239 template <class T>
240 void Vpow(int n, const Array<OneD, const T> &x, const int incx, const T f,
241  Array<OneD, T> &y, const int incy)
242 {
243  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
244  ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
245 
246  Vpow(n, &x[0], incx, f, &y[0], incy);
247 }
248 
249 /// \brief sqrt y = sqrt(x)
250 template <class T>
251 void Vsqrt(int n, const Array<OneD, const T> &x, const int incx,
252  Array<OneD, T> &y, const int incy)
253 {
254  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
255  ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
256 
257  Vsqrt(n, &x[0], incx, &y[0], incy);
258 }
259 
260 /// \brief vabs: y = |x|
261 template <class T>
262 void Vabs(int n, const Array<OneD, const T> &x, const int incx,
263  Array<OneD, T> &y, const int incy)
264 {
265  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
266  ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
267 
268  Vabs(n, &x[0], incx, &y[0], incy);
269 }
270 
271 /********** Triad routines ***********************/
272 
273 /// \brief vvtvp (vector times vector plus vector): z = w*x + y
274 template <class T>
275 void Vvtvp(int n, const Array<OneD, const T> &w, const int incw,
276  const Array<OneD, const T> &x, const int incx,
277  const Array<OneD, const T> &y, const int incy, Array<OneD, T> &z,
278  const int incz)
279 {
280  ASSERTL1(n * incw <= w.size() + w.GetOffset(), "Array out of bounds");
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 #ifdef NEKTAR_ENABLE_SIMD_VMATH
286  boost::ignore_unused(incw, incx, incy, incz);
287  ASSERTL1(incw == 1, "Simd vmath requires inc = 1");
288  ASSERTL1(incx == 1, "Simd vmath requires inc = 1");
289  ASSERTL1(incy == 1, "Simd vmath requires inc = 1");
290  ASSERTL1(incz == 1, "Simd vmath requires inc = 1");
291  SIMD::Vvtvp(n, &w[0], &x[0], &y[0], &z[0]);
292 #else
293  Vvtvp(n, &w[0], incw, &x[0], incx, &y[0], incy, &z[0], incz);
294 #endif
295 }
296 
297 template <class T>
298 void Vvtvp(int n, const Array<TwoD, NekDouble>::const_reference &w,
299  const int incw, const Array<OneD, const T> &x, const int incx,
300  const Array<OneD, const T> &y, const int incy, Array<OneD, T> &z,
301  const int incz)
302 {
303  ASSERTL1(n * incw <= w.size(), "Array out of bounds");
304  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
305  ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
306  ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
307 
308  Vvtvp(n, w.origin(), incw, &x[0], incx, &y[0], incy, &z[0], incz);
309 }
310 
311 /// \brief svtvp (scalar times vector plus vector): z = alpha*x + y
312 template <class T>
313 void Svtvp(int n, const T alpha, const Array<OneD, const T> &x, const int incx,
314  const Array<OneD, const T> &y, const int incy, Array<OneD, T> &z,
315  const int incz)
316 {
317  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
318  ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
319  ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
320 
321  Svtvp(n, alpha, &x[0], incx, &y[0], incy, &z[0], incz);
322 }
323 
324 /// \brief svtvp (scalar times vector plus vector): z = alpha*x + y
325 template <class T>
326 void Svtvm(int n, const T alpha, const Array<OneD, const T> &x, const int incx,
327  const Array<OneD, const T> &y, const int incy, Array<OneD, T> &z,
328  const int incz)
329 {
330  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
331  ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
332  ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
333 
334  Svtvm(n, alpha, &x[0], incx, &y[0], incy, &z[0], incz);
335 }
336 
337 /// \brief vvtvm (vector times vector minus vector): z = w*x - y
338 template <class T>
339 void Vvtvm(int n, const Array<OneD, const T> &w, const int incw,
340  const Array<OneD, const T> &x, const int incx,
341  const Array<OneD, const T> &y, const int incy, Array<OneD, T> &z,
342  const int incz)
343 {
344  ASSERTL1(n * incw <= w.size() + w.GetOffset(), "Array out of bounds");
345  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
346  ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
347  ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
348 
349 #ifdef NEKTAR_ENABLE_SIMD_VMATH
350  boost::ignore_unused(incw, incx, incy, incz);
351  ASSERTL1(incw == 1, "Simd vmath requires inc = 1");
352  ASSERTL1(incx == 1, "Simd vmath requires inc = 1");
353  ASSERTL1(incy == 1, "Simd vmath requires inc = 1");
354  ASSERTL1(incz == 1, "Simd vmath requires inc = 1");
355  SIMD::Vvtvm(n, &w[0], &x[0], &y[0], &z[0]);
356 #else
357  Vvtvm(n, &w[0], incw, &x[0], incx, &y[0], incy, &z[0], incz);
358 #endif
359 }
360 
361 /// \brief vvtvvtp (vector times vector plus vector times vector): z = v*w + x*y
362 template <class T>
363 void Vvtvvtp(int n, const Array<OneD, const T> &v, int incv,
364  const Array<OneD, const T> &w, int incw,
365  const Array<OneD, const T> &x, int incx,
366  const Array<OneD, const T> &y, int incy, Array<OneD, T> &z,
367  int incz)
368 {
369  ASSERTL1(n * incv <= v.size() + v.GetOffset(), "Array out of bounds");
370  ASSERTL1(n * incw <= w.size() + w.GetOffset(), "Array out of bounds");
371  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
372  ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
373  ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
374 
375 #ifdef NEKTAR_ENABLE_SIMD_VMATH
376  boost::ignore_unused(incv, incw, incx, incy, incz);
377  ASSERTL1(incw == 1, "Simd vmath requires inc = 1");
378  ASSERTL1(incx == 1, "Simd vmath requires inc = 1");
379  ASSERTL1(incy == 1, "Simd vmath requires inc = 1");
380  ASSERTL1(incz == 1, "Simd vmath requires inc = 1");
381  SIMD::Vvtvvtp(n, &v[0], &w[0], &x[0], &y[0], &z[0]);
382 #else
383  Vvtvvtp(n, &v[0], incv, &w[0], incw, &x[0], incx, &y[0], incy, &z[0], incz);
384 #endif
385 }
386 
387 /// \brief vvtvvtm (vector times vector minus vector times vector): z = v*w -
388 /// x*y
389 template <class T>
390 void Vvtvvtm(int n, const Array<OneD, const T> &v, int incv,
391  const Array<OneD, const T> &w, int incw,
392  const Array<OneD, const T> &x, int incx,
393  const Array<OneD, const T> &y, int incy, Array<OneD, T> &z,
394  int incz)
395 {
396  ASSERTL1(n * incv <= v.size() + v.GetOffset(), "Array out of bounds");
397  ASSERTL1(n * incw <= w.size() + w.GetOffset(), "Array out of bounds");
398  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
399  ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
400  ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
401 
402 #ifdef NEKTAR_ENABLE_SIMD_VMATH
403  boost::ignore_unused(incv, incw, incx, incy, incz);
404  ASSERTL1(incw == 1, "Simd vmath requires inc = 1");
405  ASSERTL1(incx == 1, "Simd vmath requires inc = 1");
406  ASSERTL1(incy == 1, "Simd vmath requires inc = 1");
407  ASSERTL1(incz == 1, "Simd vmath requires inc = 1");
408  SIMD::Vvtvvtm(n, &v[0], &w[0], &x[0], &y[0], &z[0]);
409 #else
410  Vvtvvtm(n, &v[0], incv, &w[0], incw, &x[0], incx, &y[0], incy, &z[0], incz);
411 #endif
412 }
413 
414 /// \brief svtsvtp (scalar times vector plus scalar times vector): z = alpha*x +
415 /// beta*y
416 template <class T>
417 void Svtsvtp(int n, const T alpha, const Array<OneD, const T> &x,
418  const int incx, const T beta, const Array<OneD, const T> &y,
419  const int incy, Array<OneD, T> &z, const int incz)
420 {
421  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
422  ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
423  ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
424 
425  Svtsvtp(n, alpha, &x[0], incx, beta, &y[0], incy, &z[0], incz);
426 }
427 
428 /************ Misc routine from Veclib (and extras) ************/
429 
430 /// \brief Gather vector z[i] = x[y[i]]
431 template <class T, class I,
432  typename = typename std::enable_if<std::is_floating_point<T>::value &&
433  std::is_integral<I>::value>::type>
434 void Gathr(I n, const Array<OneD, const T> &x, const Array<OneD, I> &y,
435  Array<OneD, T> &z)
436 {
437  ASSERTL1(n <= y.size() + y.GetOffset(), "Array out of bounds");
438  ASSERTL1(n <= z.size() + z.GetOffset(), "Array out of bounds");
439 
440 #ifdef NEKTAR_ENABLE_SIMD_VMATH
441  SIMD::Gathr(n, &x[0], &y[0], &z[0]);
442 #else
443  Gathr(n, &x[0], &y[0], &z[0]);
444 #endif
445 }
446 
447 /// \brief Scatter vector z[y[i]] = x[i]
448 template <class T>
449 void Scatr(int n, const Array<OneD, const T> &x,
450  const Array<OneD, const int> &y, Array<OneD, T> &z)
451 {
452  ASSERTL1(n <= x.size() + x.GetOffset(), "Array out of bounds");
453  ASSERTL1(n <= y.size() + y.GetOffset(), "Array out of bounds");
454 
455  Scatr(n, &x[0], &y[0], &z[0]);
456 }
457 
458 /// \brief Assemble z[y[i]] += x[i]; z should be zero'd first
459 template <class T>
460 void Assmb(int n, const Array<OneD, T> &x, const Array<OneD, int> &y,
461  Array<OneD, T> &z)
462 {
463  ASSERTL1(n <= x.size() + x.GetOffset(), "Array out of bounds");
464  ASSERTL1(n <= y.size() + y.GetOffset(), "Array out of bounds");
465 
466  Assmb(n, &x[0], &y[0], &z[0]);
467 }
468 
469 /************* Reduction routines *****************/
470 
471 /// \brief Subtract return sum(x)
472 template <class T> T Vsum(int n, const Array<OneD, const T> &x, const int incx)
473 {
474  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
475 
476  return Vsum(n, &x[0], incx);
477 }
478 
479 /// \brief Return the index of the maximum element in x
480 template <class T>
481 int Imax(int n, const Array<OneD, const T> &x, const int incx)
482 {
483  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
484 
485  return Imax(n, &x[0], incx);
486 }
487 
488 /// \brief Return the maximum element in x -- called vmax to avoid
489 /// conflict with max
490 template <class T> T Vmax(int n, const Array<OneD, const T> &x, const int incx)
491 {
492  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
493 
494  return Vmax(n, &x[0], incx);
495 }
496 
497 /// \brief Return the index of the maximum absolute element in x
498 template <class T>
499 int Iamax(int n, const Array<OneD, const T> &x, const int incx)
500 {
501  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
502 
503  return Iamax(n, &x[0], incx);
504 }
505 
506 /// \brief Return the maximum absolute element in x
507 /// called vamax to avoid conflict with max
508 template <class T> T Vamax(int n, const Array<OneD, const T> &x, const int incx)
509 {
510  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
511 
512  return Vamax(n, &x[0], incx);
513 }
514 
515 /// \brief Return the index of the minimum element in x
516 template <class T>
517 int Imin(int n, const Array<OneD, const T> &x, const int incx)
518 {
519  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
520 
521  return Imin(n, &x[0], incx);
522 }
523 
524 /// \brief Return the minimum element in x - called vmin to avoid
525 /// conflict with min
526 template <class T> T Vmin(int n, const Array<OneD, const T> &x, const int incx)
527 {
528  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
529 
530  return Vmin(n, &x[0], incx);
531 }
532 
533 /// \brief Return number of NaN elements of x
534 template <class T>
535 int Nnan(int n, const Array<OneD, const T> &x, const int incx)
536 {
537  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
538 
539  return Nnan(n, &x[0], incx);
540 }
541 
542 /// \brief
543 template <class T>
544 T Dot(int n, const Array<OneD, const T> &w, const Array<OneD, const T> &x)
545 {
546  ASSERTL1(n <= w.size() + w.GetOffset(), "Array out of bounds");
547  ASSERTL1(n <= x.size() + x.GetOffset(), "Array out of bounds");
548 
549  return Dot(n, &w[0], &x[0]);
550 }
551 
552 /// \brief
553 template <class T>
554 T Dot(int n, const Array<OneD, const T> &w, const int incw,
555  const Array<OneD, const T> &x, const int incx)
556 {
557  ASSERTL1(n * incw <= w.size() + w.GetOffset(), "Array out of bounds");
558  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
559 
560  return Dot(n, &w[0], incw, &x[0], incx);
561 }
562 
563 /// \brief
564 template <class T>
565 T Dot2(int n, const Array<OneD, const T> &w, const Array<OneD, const T> &x,
566  const Array<OneD, const int> &y)
567 {
568  ASSERTL1(n <= w.size() + w.GetOffset(), "Array out of bounds");
569  ASSERTL1(n <= x.size() + x.GetOffset(), "Array out of bounds");
570  ASSERTL1(n <= y.size() + y.GetOffset(), "Array out of bounds");
571 
572  return Dot2(n, &w[0], &x[0], &y[0]);
573 }
574 
575 /// \brief
576 template <class T>
577 T Ddot(int n, const Array<OneD, const T> &w, const int incw,
578  const Array<OneD, const T> &x, const int incx,
579  const Array<OneD, const int> &y, const int incy)
580 {
581  ASSERTL1(n * incw <= w.size() + w.GetOffset(), "Array out of bounds");
582  ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
583  ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
584 
585  return Dot2(n, &w[0], incw, &x[0], incx, &y[0], incy);
586 }
587 
588 /********** Memory routines ***********************/
589 
590 template <class T>
591 void Vcopy(int n, const Array<OneD, const T> &x, int incx, Array<OneD, T> &y,
592  int const incy)
593 {
594  ASSERTL1(static_cast<unsigned int>(std::abs(n * incx)) <=
595  x.size() + x.GetOffset(),
596  "Array out of bounds");
597  ASSERTL1(static_cast<unsigned int>(std::abs(n * incy)) <=
598  y.size() + y.GetOffset(),
599  "Array out of bounds");
600 
601  Vcopy(n, &x[0], incx, &y[0], incy);
602 }
603 
604 template <class T>
605 void Reverse(int n, const Array<OneD, const T> &x, int incx, Array<OneD, T> &y,
606  int const incy)
607 {
608  ASSERTL1(static_cast<unsigned int>(std::abs(n * incx)) <=
609  x.size() + x.GetOffset(),
610  "Array out of bounds");
611  ASSERTL1(static_cast<unsigned int>(std::abs(n * incy)) <=
612  y.size() + y.GetOffset(),
613  "Array out of bounds");
614 
615  Reverse(n, &x[0], incx, &y[0], incy);
616 }
617 
618 } // namespace Vmath
619 #endif // VECTORMATHARRAY_HPP
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:249
Nektar::LibUtilities::beta
@ beta
Gauss Radau pinned at x=-1,.
Definition: PointsType.h:61
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:2
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:263
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:51
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:311
Vmath::SIMD::Vvtvvtm
void Vvtvvtm(const size_t n, const T *v, const T *w, const T *x, const T *y, T *z)
vvtvvtm (vector times vector minus vector times vector):
Definition: VmathSIMD.hpp:418
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:476
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:360
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:157
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:534
Vmath::Ssub
void Ssub(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Substract vector y = alpha - x.
Definition: Vmath.cpp:458
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)
svtvvtp (scalar times vector plus scalar times vector):
Definition: Vmath.cpp:751
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:209
Vmath::Vlog
void Vlog(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:114
Vmath::Vexp
void Vexp(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:125
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:622
Vmath::Vabs
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
Definition: Vmath.cpp:553
Vmath::Dot2
T Dot2(int n, const T *w, const T *x, const int *y)
dot2 (vector times vector times vector): z = w*x*y
Definition: Vmath.cpp:1147
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:518
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:1050
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:577
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:574
Vmath::Vsum
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:895
Vmath::Scatr
void Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
Definition: Vmath.cpp:822
Vmath::Dot
T Dot(int n, const T *w, const T *x)
dot (vector times vector): z = w*x
Definition: Vmath.cpp:1100
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:805
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:862
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 minus vector): z = alpha*x - y
Definition: Vmath.cpp:664
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:359
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 minus vector): z = w*x - y
Definition: Vmath.cpp:598
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:721
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:248
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:324
Vmath::Imax
int Imax(int n, const T *x, const int incx)
Return the index of the maximum element in x.
Definition: Vmath.cpp:918
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:284
Vmath::Imin
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
Definition: Vmath.cpp:1023
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492
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:155
Vmath::Nnan
int Nnan(int n, const T *x, const int incx)
Return number of NaN elements of x.
Definition: Vmath.cpp:1076
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:999
Vmath::Sadd
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add scalar y = alpha + x.
Definition: Vmath.cpp:384
Vmath::Reverse
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1286
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:945
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:692
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255
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:419
Vmath::Vpow
void Vpow(int n, const T *x, const int incx, const T f, T *y, const int incy)
Definition: Vmath.hpp:136
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:971
tinysimd::abs
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:298
Generated on Sat Mar 25 2023 12:35:53 for Nektar++ by doxygen 1.9.1