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 arguments
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
42namespace Vmath
43{
44using namespace Nektar;
45
46/***************** Math routines ***************/
47/// \brief Fill a vector with a constant value
48template <class T>
49void 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
57template <class T>
58void 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
67template <class T>
68void Vmul(int n, const Array<OneD, const T> &x, [[maybe_unused]] const int incx,
69 const Array<OneD, const T> &y, [[maybe_unused]] const int incy,
70 Array<OneD, T> &z, [[maybe_unused]] 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 ASSERTL1(incx == 1, "Simd vmath requires inc = 1");
78 ASSERTL1(incy == 1, "Simd vmath requires inc = 1");
79 ASSERTL1(incz == 1, "Simd vmath requires inc = 1");
80 SIMD::Vmul(n, &x[0], &y[0], &z[0]);
81#else
82 Vmul(n, &x[0], incx, &y[0], incy, &z[0], incz);
83#endif
84}
85
86template <class T>
87void Vmul(int n, const Array<TwoD, NekDouble>::const_reference &x,
88 const int incx, const Array<OneD, const T> &y, const int incy,
89 Array<OneD, T> &z, const int incz)
90{
91 ASSERTL1(n * incx <= x.size(), "Array out of bounds");
92 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
93 ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
94
95 Vmul(n, x.origin(), incx, &y[0], incy, &z[0], incz);
96}
97
98/// \brief Scalar multiply y = alpha*x
99
100template <class T>
101void Smul(int n, const T alpha, const Array<OneD, const T> &x, const int incx,
102 Array<OneD, T> &y, const int incy)
103{
104 ASSERTL1(static_cast<unsigned int>(n * incx) <= x.size() + x.GetOffset(),
105 "Array out of bounds");
106 ASSERTL1(static_cast<unsigned int>(n * incy) <= y.size() + y.GetOffset(),
107 "Array out of bounds");
108
109 Smul(n, alpha, &x[0], incx, &y[0], incy);
110}
111
112/// \brief Multiply vector z = x/y
113template <class T>
114void Vdiv(int n, const Array<OneD, const T> &x, const int incx,
115 const Array<OneD, const T> &y, const int incy, Array<OneD, T> &z,
116 const int incz)
117{
118 ASSERTL1(static_cast<unsigned int>(n * incx) <= x.size() + x.GetOffset(),
119 "Array out of bounds");
120 ASSERTL1(static_cast<unsigned int>(n * incy) <= y.size() + y.GetOffset(),
121 "Array out of bounds");
122 ASSERTL1(static_cast<unsigned int>(n * incz) <= z.size() + z.GetOffset(),
123 "Array out of bounds");
124
125 Vdiv(n, &x[0], incx, &y[0], incy, &z[0], incz);
126}
127
128/// \brief Scalar multiply y = alpha/x
129template <class T>
130void Sdiv(int n, const T alpha, const Array<OneD, const T> &x, const int incx,
131 Array<OneD, T> &y, const int incy)
132{
133 ASSERTL1(static_cast<unsigned int>(n * incx) <= x.size() + x.GetOffset(),
134 "Array out of bounds");
135 ASSERTL1(static_cast<unsigned int>(n * incy) <= y.size() + y.GetOffset(),
136 "Array out of bounds");
137
138 Sdiv(n, alpha, &x[0], incx, &y[0], incy);
139}
140
141/// \brief Add vector z = x+y
142template <class T>
143void Vadd(int n, const Array<OneD, const T> &x, [[maybe_unused]] const int incx,
144 const Array<OneD, const T> &y, [[maybe_unused]] const int incy,
145 Array<OneD, T> &z, [[maybe_unused]] const int incz)
146{
147 ASSERTL1(static_cast<unsigned int>(n * incx) <= x.size() + x.GetOffset(),
148 "Array out of bounds");
149 ASSERTL1(static_cast<unsigned int>(n * incy) <= y.size() + y.GetOffset(),
150 "Array out of bounds");
151 ASSERTL1(static_cast<unsigned int>(n * incz) <= z.size() + z.GetOffset(),
152 "Array out of bounds");
153
154#ifdef NEKTAR_ENABLE_SIMD_VMATH
155 ASSERTL1(incx == 1, "Simd vmath requires inc = 1");
156 ASSERTL1(incy == 1, "Simd vmath requires inc = 1");
157 ASSERTL1(incz == 1, "Simd vmath requires inc = 1");
158 SIMD::Vadd(n, &x[0], &y[0], &z[0]);
159#else
160 Vadd(n, &x[0], incx, &y[0], incy, &z[0], incz);
161#endif
162}
163
164/// \brief Add vector y = alpha + x
165template <class T>
166void Sadd(int n, const T alpha, const Array<OneD, const T> &x, const int incx,
167 Array<OneD, T> &y, const int incy)
168{
169
170 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
171 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
172
173 Sadd(n, alpha, &x[0], incx, &y[0], incy);
174}
175
176/// \brief Subtract vector z = x-y
177template <class T>
178void Vsub(int n, const Array<OneD, const T> &x, const int incx,
179 const Array<OneD, const T> &y, const int incy, Array<OneD, T> &z,
180 const int incz)
181{
182 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
183 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
184 ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
185
186 Vsub(n, &x[0], incx, &y[0], incy, &z[0], incz);
187}
188
189/// \brief Add vector y = alpha - x
190template <class T>
191void Ssub(int n, const T alpha, const Array<OneD, const T> &x, const int incx,
192 Array<OneD, T> &y, const int incy)
193{
194
195 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
196 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
197
198 Ssub(n, alpha, &x[0], incx, &y[0], incy);
199}
200
201/// \brief Zero vector
202template <class T> void Zero(int n, Array<OneD, T> &x, const int incx)
203{
204 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
205
206 Zero(n, &x[0], incx);
207}
208
209/// \brief Negate x = -x
210template <class T> void Neg(int n, Array<OneD, T> &x, const int incx)
211{
212 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
213
214 Neg(n, &x[0], incx);
215}
216
217/// \brief log y = log(x)
218template <class T>
219void Vlog(int n, const Array<OneD, const T> &x, const int incx,
220 Array<OneD, T> &y, const int incy)
221{
222 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
223 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
224
225 Vlog(n, &x[0], incx, &y[0], incy);
226}
227
228/// \brief exp y = exp(x)
229template <class T>
230void 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/// \brief pow y = pow(x, f)
240template <class T>
241void Vpow(int n, const Array<OneD, const T> &x, const int incx, const T f,
242 Array<OneD, T> &y, const int incy)
243{
244 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
245 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
246
247 Vpow(n, &x[0], incx, f, &y[0], incy);
248}
249
250/// \brief sqrt y = sqrt(x)
251template <class T>
252void Vsqrt(int n, const Array<OneD, const T> &x, const int incx,
253 Array<OneD, T> &y, const int incy)
254{
255 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
256 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
257
258 Vsqrt(n, &x[0], incx, &y[0], incy);
259}
260
261/// \brief vabs: y = |x|
262template <class T>
263void Vabs(int n, const Array<OneD, const T> &x, const int incx,
264 Array<OneD, T> &y, const int incy)
265{
266 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
267 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
268
269 Vabs(n, &x[0], incx, &y[0], incy);
270}
271
272/********** Triad routines ***********************/
273
274/// \brief vvtvp (vector times vector plus vector): z = w*x + y
275template <class T>
276void Vvtvp(int n, const Array<OneD, const T> &w,
277 [[maybe_unused]] const int incw, const Array<OneD, const T> &x,
278 [[maybe_unused]] const int incx, const Array<OneD, const T> &y,
279 [[maybe_unused]] const int incy, Array<OneD, T> &z,
280 [[maybe_unused]] const int incz)
281{
282 ASSERTL1(n * incw <= w.size() + w.GetOffset(), "Array out of bounds");
283 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
284 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
285 ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
286
287#ifdef NEKTAR_ENABLE_SIMD_VMATH
288 ASSERTL1(incw == 1, "Simd vmath requires inc = 1");
289 ASSERTL1(incx == 1, "Simd vmath requires inc = 1");
290 ASSERTL1(incy == 1, "Simd vmath requires inc = 1");
291 ASSERTL1(incz == 1, "Simd vmath requires inc = 1");
292 SIMD::Vvtvp(n, &w[0], &x[0], &y[0], &z[0]);
293#else
294 Vvtvp(n, &w[0], incw, &x[0], incx, &y[0], incy, &z[0], incz);
295#endif
296}
297
298/// \brief Vvtvp (vector times vector plus vector): z = w*x + y
299template <class T>
300void Vvtvp(int n, const Array<TwoD, NekDouble>::const_reference &w,
301 const int incw, const Array<OneD, const T> &x, const int incx,
302 const Array<OneD, const T> &y, const int incy, Array<OneD, T> &z,
303 const int incz)
304{
305 ASSERTL1(n * incw <= w.size(), "Array out of bounds");
306 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
307 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
308 ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
309
310 Vvtvp(n, w.origin(), incw, &x[0], incx, &y[0], incy, &z[0], incz);
311}
312
313/// \brief vvtvm (vector times vector minus vector): z = w*x - y
314template <class T>
315void Vvtvm(int n, const Array<OneD, const T> &w,
316 [[maybe_unused]] const int incw, const Array<OneD, const T> &x,
317 [[maybe_unused]] const int incx, const Array<OneD, const T> &y,
318 [[maybe_unused]] const int incy, Array<OneD, T> &z,
319 [[maybe_unused]] const int incz)
320{
321 ASSERTL1(n * incw <= w.size() + w.GetOffset(), "Array out of bounds");
322 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
323 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
324 ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
325
326#ifdef NEKTAR_ENABLE_SIMD_VMATH
327 ASSERTL1(incw == 1, "Simd vmath requires inc = 1");
328 ASSERTL1(incx == 1, "Simd vmath requires inc = 1");
329 ASSERTL1(incy == 1, "Simd vmath requires inc = 1");
330 ASSERTL1(incz == 1, "Simd vmath requires inc = 1");
331 SIMD::Vvtvm(n, &w[0], &x[0], &y[0], &z[0]);
332#else
333 Vvtvm(n, &w[0], incw, &x[0], incx, &y[0], incy, &z[0], incz);
334#endif
335}
336
337/// \brief svtvp (scalar times vector plus vector): z = alpha*x + y
338template <class T>
339void Svtvp(int n, const T alpha, const Array<OneD, const T> &x, const int incx,
340 const Array<OneD, const T> &y, const int incy, Array<OneD, T> &z,
341 const int incz)
342{
343 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
344 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
345 ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
346
347 Svtvp(n, alpha, &x[0], incx, &y[0], incy, &z[0], incz);
348}
349
350/// \brief svtvm (scalar times vector minus vector): z = alpha*x - y
351template <class T>
352void Svtvm(int n, const T alpha, const Array<OneD, const T> &x, const int incx,
353 const Array<OneD, const T> &y, const int incy, Array<OneD, T> &z,
354 const int incz)
355{
356 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
357 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
358 ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
359
360 Svtvm(n, alpha, &x[0], incx, &y[0], incy, &z[0], incz);
361}
362
363/// \brief vvtvvtp (vector times vector plus vector times vector): z = v*w + x*y
364template <class T>
365void Vvtvvtp(int n, const Array<OneD, const T> &v, int incv,
366 const Array<OneD, const T> &w, [[maybe_unused]] int incw,
367 const Array<OneD, const T> &x, [[maybe_unused]] int incx,
368 const Array<OneD, const T> &y, [[maybe_unused]] int incy,
369 Array<OneD, T> &z, [[maybe_unused]] int incz)
370{
371 ASSERTL1(n * incv <= v.size() + v.GetOffset(), "Array out of bounds");
372 ASSERTL1(n * incw <= w.size() + w.GetOffset(), "Array out of bounds");
373 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
374 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
375 ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
376
377#ifdef NEKTAR_ENABLE_SIMD_VMATH
378 ASSERTL1(incw == 1, "Simd vmath requires inc = 1");
379 ASSERTL1(incx == 1, "Simd vmath requires inc = 1");
380 ASSERTL1(incy == 1, "Simd vmath requires inc = 1");
381 ASSERTL1(incz == 1, "Simd vmath requires inc = 1");
382 SIMD::Vvtvvtp(n, &v[0], &w[0], &x[0], &y[0], &z[0]);
383#else
384 Vvtvvtp(n, &v[0], incv, &w[0], incw, &x[0], incx, &y[0], incy, &z[0], incz);
385#endif
386}
387
388/// \brief vvtvvtm (vector times vector minus vector times vector): z = v*w -
389/// x*y
390template <class T>
391void Vvtvvtm(int n, const Array<OneD, const T> &v, [[maybe_unused]] int incv,
392 const Array<OneD, const T> &w, [[maybe_unused]] int incw,
393 const Array<OneD, const T> &x, [[maybe_unused]] int incx,
394 const Array<OneD, const T> &y, [[maybe_unused]] int incy,
395 Array<OneD, T> &z, [[maybe_unused]] int incz)
396{
397 ASSERTL1(n * incv <= v.size() + v.GetOffset(), "Array out of bounds");
398 ASSERTL1(n * incw <= w.size() + w.GetOffset(), "Array out of bounds");
399 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
400 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
401 ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
402
403#ifdef NEKTAR_ENABLE_SIMD_VMATH
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
416template <class T>
417void 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]]
431template <class T, class I,
432 typename = typename std::enable_if<std::is_floating_point<T>::value &&
433 std::is_integral<I>::value>::type>
434void 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]
448template <class T>
449void 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
459template <class T>
460void 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)
472template <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
480template <class T>
481int 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
490template <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
498template <class T>
499int 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
508template <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
516template <class T>
517int 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
526template <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
534template <class T>
535int 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 dot product
543template <class T>
544T 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 dot product
553template <class T>
554T 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 dot product
564template <class T>
565T 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 dot product
576template <class T>
577T 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// \brief copy one vector to another
591template <class T>
592void Vcopy(int n, const Array<OneD, const T> &x, int incx, Array<OneD, T> &y,
593 int const incy)
594{
595 ASSERTL1(static_cast<unsigned int>(std::abs(n * incx)) <=
596 x.size() + x.GetOffset(),
597 "Array out of bounds");
598 ASSERTL1(static_cast<unsigned int>(std::abs(n * incy)) <=
599 y.size() + y.GetOffset(),
600 "Array out of bounds");
601
602 Vcopy(n, &x[0], incx, &y[0], incy);
603}
604
605// \brief reverse the ordering of vector to another
606template <class T>
607void Reverse(int n, const Array<OneD, const T> &x, int incx, Array<OneD, T> &y,
608 int const incy)
609{
610 ASSERTL1(static_cast<unsigned int>(std::abs(n * incx)) <=
611 x.size() + x.GetOffset(),
612 "Array out of bounds");
613 ASSERTL1(static_cast<unsigned int>(std::abs(n * incy)) <=
614 y.size() + y.GetOffset(),
615 "Array out of bounds");
616
617 Reverse(n, &x[0], incx, &y[0], incy);
618}
619
620} // namespace Vmath
621#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:242
Nektar::LibUtilities::beta
@ beta
Gauss Radau pinned at x=-1,.
Definition: PointsType.h:59
Nektar::UnitTests::w
std::vector< double > w(NPUPPER)
Nektar::UnitTests::z
std::vector< double > z(NPUPPER)
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:259
Vmath::SIMD::Vadd
void Vadd(const size_t n, const T *x, const T *y, T *z)
Add vector z = x + y.
Definition: VmathSIMD.hpp:47
Vmath::SIMD::Vvtvm
void Vvtvm(const size_t n, const T *w, const T *x, const T *y, T *z)
vvtvm (vector times vector minus vector): z = w*x - y
Definition: VmathSIMD.hpp:307
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:414
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:472
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:356
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:153
Vmath
Definition: Vmath.cpp:38
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.hpp:340
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.hpp:248
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)
Svtsvtp (scalar times vector plus scalar times vector):
Definition: Vmath.hpp:473
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.hpp:72
Vmath::Gathr
void Gathr(I n, const T *x, const I *y, T *z)
Gather vector z[i] = x[y[i]].
Definition: Vmath.hpp:507
Vmath::Vlog
void Vlog(int n, const T *x, const int incx, T *y, const int incy)
log y = log(x)
Definition: Vmath.hpp:303
Vmath::Vexp
void Vexp(int n, const T *x, const int incx, T *y, const int incy)
exp y = exp(x)
Definition: Vmath.hpp:315
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.hpp:396
Vmath::Vabs
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
Definition: Vmath.hpp:352
Vmath::Dot2
T Dot2(int n, const T *w, const T *x, const int *y)
dot product
Definition: Vmath.hpp:790
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.hpp:292
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.hpp:725
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)
dot product
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.hpp:366
Vmath::Vsum
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.hpp:608
Vmath::Scatr
void Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
Definition: Vmath.hpp:539
Vmath::Dot
T Dot(int n, const T *w, const T *x)
dot product
Definition: Vmath.hpp:761
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.hpp:577
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)
Svtvm (scalar times vector minus vector): z = alpha*x - y.
Definition: Vmath.hpp:424
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.hpp:180
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.hpp:381
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.hpp:456
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.hpp:100
Vmath::Sdiv
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/x.
Definition: Vmath.hpp:154
Vmath::Imax
int Imax(int n, const T *x, const int incx)
Return the index of the maximum element in x.
Definition: Vmath.hpp:623
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.hpp:126
Vmath::Imin
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
Definition: Vmath.hpp:704
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.hpp:54
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:154
Vmath::Nnan
int Nnan(int n, const T *x, const int incx)
Return number of NaN elements of x.
Definition: Vmath.hpp:743
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.hpp:685
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.hpp:194
Vmath::Reverse
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:844
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.hpp:644
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.hpp:439
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
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.hpp:220
Vmath::Vpow
void Vpow(int n, const T *x, const int incx, const T f, T *y, const int incy)
pow y = pow(x, f)
Definition: Vmath.hpp:327
Vmath::Iamax
int Iamax(int n, const T *x, const int incx)
Return the index of the maximum absolute element in x.
Definition: Vmath.hpp:662
tinysimd::abs
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:298
Generated on Wed Jan 31 2024 09:35:54 for Nektar++ by doxygen 1.9.4