Nektar++
Loading...
Searching...
No Matches
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/SimdLib/tinysimd.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 Vmul(n, &x[0], incx, &y[0], incy, &z[0], incz);
77}
78
79template <class T>
80void Vmul(int n, const typename Array<TwoD, T>::const_reference &x,
81 const int incx, const Array<OneD, const T> &y, const int incy,
82 Array<OneD, T> &z, const int incz)
83{
84 ASSERTL1(n * incx <= x.size(), "Array out of bounds");
85 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
86 ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
87
88 Vmul(n, x.origin(), incx, &y[0], incy, &z[0], incz);
89}
90
91/// \brief Scalar multiply y = alpha*x
92
93template <class T>
94void Smul(int n, const T alpha, const Array<OneD, const T> &x, const int incx,
95 Array<OneD, T> &y, const int incy)
96{
97 ASSERTL1(static_cast<unsigned int>(n * incx) <= x.size() + x.GetOffset(),
98 "Array out of bounds");
99 ASSERTL1(static_cast<unsigned int>(n * incy) <= y.size() + y.GetOffset(),
100 "Array out of bounds");
101
102 Smul(n, alpha, &x[0], incx, &y[0], incy);
103}
104
105/// \brief Multiply vector z = x/y
106template <class T>
107void Vdiv(int n, const Array<OneD, const T> &x, const int incx,
108 const Array<OneD, const T> &y, const int incy, Array<OneD, T> &z,
109 const int incz)
110{
111 ASSERTL1(static_cast<unsigned int>(n * incx) <= x.size() + x.GetOffset(),
112 "Array out of bounds");
113 ASSERTL1(static_cast<unsigned int>(n * incy) <= y.size() + y.GetOffset(),
114 "Array out of bounds");
115 ASSERTL1(static_cast<unsigned int>(n * incz) <= z.size() + z.GetOffset(),
116 "Array out of bounds");
117
118 Vdiv(n, &x[0], incx, &y[0], incy, &z[0], incz);
119}
120
121/// \brief Scalar multiply y = alpha/x
122template <class T>
123void Sdiv(int n, const T alpha, const Array<OneD, const T> &x, const int incx,
124 Array<OneD, T> &y, const int incy)
125{
126 ASSERTL1(static_cast<unsigned int>(n * incx) <= x.size() + x.GetOffset(),
127 "Array out of bounds");
128 ASSERTL1(static_cast<unsigned int>(n * incy) <= y.size() + y.GetOffset(),
129 "Array out of bounds");
130
131 Sdiv(n, alpha, &x[0], incx, &y[0], incy);
132}
133
134/// \brief Add vector z = x+y
135template <class T>
136void Vadd(int n, const Array<OneD, const T> &x, [[maybe_unused]] const int incx,
137 const Array<OneD, const T> &y, [[maybe_unused]] const int incy,
138 Array<OneD, T> &z, [[maybe_unused]] const int incz)
139{
140 ASSERTL1(static_cast<unsigned int>(n * incx) <= x.size() + x.GetOffset(),
141 "Array out of bounds");
142 ASSERTL1(static_cast<unsigned int>(n * incy) <= y.size() + y.GetOffset(),
143 "Array out of bounds");
144 ASSERTL1(static_cast<unsigned int>(n * incz) <= z.size() + z.GetOffset(),
145 "Array out of bounds");
146
147 Vadd(n, &x[0], incx, &y[0], incy, &z[0], incz);
148}
149
150/// \brief Add vector y = alpha + x
151template <class T>
152void Sadd(int n, const T alpha, const Array<OneD, const T> &x, const int incx,
153 Array<OneD, T> &y, const int incy)
154{
155
156 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
157 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
158
159 Sadd(n, alpha, &x[0], incx, &y[0], incy);
160}
161
162/// \brief Subtract vector z = x-y
163template <class T>
164void Vsub(int n, const Array<OneD, const T> &x, const int incx,
165 const Array<OneD, const T> &y, const int incy, Array<OneD, T> &z,
166 const int incz)
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 ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
171
172 Vsub(n, &x[0], incx, &y[0], incy, &z[0], incz);
173}
174
175/// \brief Add vector y = alpha - x
176template <class T>
177void Ssub(int n, const T alpha, const Array<OneD, const T> &x, const int incx,
178 Array<OneD, T> &y, const int incy)
179{
180
181 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
182 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
183
184 Ssub(n, alpha, &x[0], incx, &y[0], incy);
185}
186
187/// \brief Zero vector
188template <class T> void Zero(int n, Array<OneD, T> &x, const int incx)
189{
190 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
191
192 Zero(n, &x[0], incx);
193}
194
195/// \brief Negate x = -x
196template <class T> void Neg(int n, Array<OneD, T> &x, const int incx)
197{
198 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
199
200 Neg(n, &x[0], incx);
201}
202
203/// \brief log y = log(x)
204template <class T>
205void Vlog(int n, const Array<OneD, const T> &x, const int incx,
206 Array<OneD, T> &y, const int incy)
207{
208 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
209 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
210
211 Vlog(n, &x[0], incx, &y[0], incy);
212}
213
214/// \brief exp y = exp(x)
215template <class T>
216void Vexp(int n, const Array<OneD, const T> &x, const int incx,
217 Array<OneD, T> &y, const int incy)
218{
219 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
220 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
221
222 Vexp(n, &x[0], incx, &y[0], incy);
223}
224
225/// \brief pow y = pow(x, f)
226template <class T>
227void Vpow(int n, const Array<OneD, const T> &x, const int incx, const T f,
228 Array<OneD, T> &y, const int incy)
229{
230 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
231 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
232
233 Vpow(n, &x[0], incx, f, &y[0], incy);
234}
235
236/// \brief sqrt y = sqrt(x)
237template <class T>
238void Vsqrt(int n, const Array<OneD, const T> &x, const int incx,
239 Array<OneD, T> &y, const int incy)
240{
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
244 Vsqrt(n, &x[0], incx, &y[0], incy);
245}
246
247/// \brief vabs: y = |x|
248template <class T>
249void Vabs(int n, const Array<OneD, const T> &x, const int incx,
250 Array<OneD, T> &y, const int incy)
251{
252 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
253 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
254
255 Vabs(n, &x[0], incx, &y[0], incy);
256}
257
258/********** Triad routines ***********************/
259
260/// \brief vvtvp (vector times vector plus vector): z = w*x + y
261template <class T>
262void Vvtvp(int n, const Array<OneD, const T> &w,
263 [[maybe_unused]] const int incw, const Array<OneD, const T> &x,
264 [[maybe_unused]] const int incx, const Array<OneD, const T> &y,
265 [[maybe_unused]] const int incy, Array<OneD, T> &z,
266 [[maybe_unused]] const int incz)
267{
268 ASSERTL1(n * incw <= w.size() + w.GetOffset(), "Array out of bounds");
269 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
270 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
271 ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
272
273 Vvtvp(n, &w[0], incw, &x[0], incx, &y[0], incy, &z[0], incz);
274}
275
276/// \brief Vvtvp (vector times vector plus vector): z = w*x + y
277template <class T>
278void Vvtvp(int n, const typename Array<TwoD, T>::const_reference &w,
279 const int incw, const Array<OneD, const T> &x, const int incx,
280 const Array<OneD, const T> &y, const int incy, Array<OneD, T> &z,
281 const int incz)
282{
283 ASSERTL1(n * incw <= w.size(), "Array out of bounds");
284 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
285 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
286 ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
287
288 Vvtvp(n, w.origin(), incw, &x[0], incx, &y[0], incy, &z[0], incz);
289}
290
291/// \brief vvtvm (vector times vector minus vector): z = w*x - y
292template <class T>
293void Vvtvm(int n, const Array<OneD, const T> &w,
294 [[maybe_unused]] const int incw, const Array<OneD, const T> &x,
295 [[maybe_unused]] const int incx, const Array<OneD, const T> &y,
296 [[maybe_unused]] const int incy, Array<OneD, T> &z,
297 [[maybe_unused]] const int incz)
298{
299 ASSERTL1(n * incw <= w.size() + w.GetOffset(), "Array out of bounds");
300 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
301 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
302 ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
303
304 Vvtvm(n, &w[0], incw, &x[0], incx, &y[0], incy, &z[0], incz);
305}
306
307/// \brief svtvp (scalar times vector plus vector): z = alpha*x + y
308template <class T>
309void Svtvp(int n, const T alpha, const Array<OneD, const T> &x, const int incx,
310 const Array<OneD, const T> &y, const int incy, Array<OneD, T> &z,
311 const int incz)
312{
313 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
314 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
315 ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
316
317 Svtvp(n, alpha, &x[0], incx, &y[0], incy, &z[0], incz);
318}
319
320/// \brief svtvm (scalar times vector minus vector): z = alpha*x - y
321template <class T>
322void Svtvm(int n, const T alpha, const Array<OneD, const T> &x, const int incx,
323 const Array<OneD, const T> &y, const int incy, Array<OneD, T> &z,
324 const int incz)
325{
326 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
327 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
328 ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
329
330 Svtvm(n, alpha, &x[0], incx, &y[0], incy, &z[0], incz);
331}
332
333/// \brief vvtvvtp (vector times vector plus vector times vector): z = v*w + x*y
334template <class T>
335void Vvtvvtp(int n, const Array<OneD, const T> &v, int incv,
336 const Array<OneD, const T> &w, [[maybe_unused]] int incw,
337 const Array<OneD, const T> &x, [[maybe_unused]] int incx,
338 const Array<OneD, const T> &y, [[maybe_unused]] int incy,
339 Array<OneD, T> &z, [[maybe_unused]] int incz)
340{
341 ASSERTL1(n * incv <= v.size() + v.GetOffset(), "Array out of bounds");
342 ASSERTL1(n * incw <= w.size() + w.GetOffset(), "Array out of bounds");
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 Vvtvvtp(n, &v[0], incv, &w[0], incw, &x[0], incx, &y[0], incy, &z[0], incz);
348}
349
350/// \brief vvtvvtm (vector times vector minus vector times vector): z = v*w -
351/// x*y
352template <class T>
353void Vvtvvtm(int n, const Array<OneD, const T> &v, [[maybe_unused]] int incv,
354 const Array<OneD, const T> &w, [[maybe_unused]] int incw,
355 const Array<OneD, const T> &x, [[maybe_unused]] int incx,
356 const Array<OneD, const T> &y, [[maybe_unused]] int incy,
357 Array<OneD, T> &z, [[maybe_unused]] int incz)
358{
359 ASSERTL1(n * incv <= v.size() + v.GetOffset(), "Array out of bounds");
360 ASSERTL1(n * incw <= w.size() + w.GetOffset(), "Array out of bounds");
361 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
362 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
363 ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
364
365 Vvtvvtm(n, &v[0], incv, &w[0], incw, &x[0], incx, &y[0], incy, &z[0], incz);
366}
367
368/// \brief svtsvtp (scalar times vector plus scalar times vector): z = alpha*x +
369/// beta*y
370template <class T>
371void Svtsvtp(int n, const T alpha, const Array<OneD, const T> &x,
372 const int incx, const T beta, const Array<OneD, const T> &y,
373 const int incy, Array<OneD, T> &z, const int incz)
374{
375 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
376 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
377 ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
378
379 Svtsvtp(n, alpha, &x[0], incx, beta, &y[0], incy, &z[0], incz);
380}
381
382/************ Misc routine from Veclib (and extras) ************/
383
384/// \brief Gather vector z[i] = x[y[i]]
385template <class T, class I,
386 typename = typename std::enable_if<std::is_floating_point_v<T> &&
387 std::is_integral_v<I>>::type>
388void Gathr(I n, const Array<OneD, const T> &x, const Array<OneD, I> &y,
389 Array<OneD, T> &z)
390{
391 ASSERTL1(n <= y.size() + y.GetOffset(), "Array out of bounds");
392 ASSERTL1(n <= z.size() + z.GetOffset(), "Array out of bounds");
393
394 Gathr(n, &x[0], &y[0], &z[0]);
395}
396
397/// \brief Scatter vector z[y[i]] = x[i]
398template <class T>
399void Scatr(int n, const Array<OneD, const T> &x,
400 const Array<OneD, const int> &y, Array<OneD, T> &z)
401{
402 ASSERTL1(n <= x.size() + x.GetOffset(), "Array out of bounds");
403 ASSERTL1(n <= y.size() + y.GetOffset(), "Array out of bounds");
404
405 Scatr(n, &x[0], &y[0], &z[0]);
406}
407
408/// \brief Assemble z[y[i]] += x[i]; z should be zero'd first
409template <class T>
410void Assmb(int n, const Array<OneD, T> &x, const Array<OneD, int> &y,
411 Array<OneD, T> &z)
412{
413 ASSERTL1(n <= x.size() + x.GetOffset(), "Array out of bounds");
414 ASSERTL1(n <= y.size() + y.GetOffset(), "Array out of bounds");
415
416 Assmb(n, &x[0], &y[0], &z[0]);
417}
418
419/************* Reduction routines *****************/
420
421/// \brief Subtract return sum(x)
422template <class T> T Vsum(int n, const Array<OneD, const T> &x, const int incx)
423{
424 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
425
426 return Vsum(n, &x[0], incx);
427}
428
429/// \brief Return the index of the maximum element in x
430template <class T>
431int Imax(int n, const Array<OneD, const T> &x, const int incx)
432{
433 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
434
435 return Imax(n, &x[0], incx);
436}
437
438/// \brief Return the maximum element in x -- called vmax to avoid
439/// conflict with max
440template <class T> T Vmax(int n, const Array<OneD, const T> &x, const int incx)
441{
442 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
443
444 return Vmax(n, &x[0], incx);
445}
446
447/// \brief Return the index of the maximum absolute element in x
448template <class T>
449int Iamax(int n, const Array<OneD, const T> &x, const int incx)
450{
451 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
452
453 return Iamax(n, &x[0], incx);
454}
455
456/// \brief Return the maximum absolute element in x
457/// called vamax to avoid conflict with max
458template <class T> T Vamax(int n, const Array<OneD, const T> &x, const int incx)
459{
460 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
461
462 return Vamax(n, &x[0], incx);
463}
464
465/// \brief Return the index of the minimum element in x
466template <class T>
467int Imin(int n, const Array<OneD, const T> &x, const int incx)
468{
469 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
470
471 return Imin(n, &x[0], incx);
472}
473
474/// \brief Return the minimum element in x - called vmin to avoid
475/// conflict with min
476template <class T> T Vmin(int n, const Array<OneD, const T> &x, const int incx)
477{
478 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
479
480 return Vmin(n, &x[0], incx);
481}
482
483/// \brief Return number of NaN elements of x
484template <class T>
485int Nnan(int n, const Array<OneD, const T> &x, const int incx)
486{
487 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
488
489 return Nnan(n, &x[0], incx);
490}
491
492/// \brief dot product
493template <class T>
494T Dot(int n, const Array<OneD, const T> &w, const Array<OneD, const T> &x)
495{
496 ASSERTL1(n <= w.size() + w.GetOffset(), "Array out of bounds");
497 ASSERTL1(n <= x.size() + x.GetOffset(), "Array out of bounds");
498
499 return Dot(n, &w[0], &x[0]);
500}
501
502/// \brief dot product
503template <class T>
504T Dot(int n, const Array<OneD, const T> &w, const int incw,
505 const Array<OneD, const T> &x, const int incx)
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
510 return Dot(n, &w[0], incw, &x[0], incx);
511}
512
513/// \brief dot product
514template <class T>
515T Dot2(int n, const Array<OneD, const T> &w, const Array<OneD, const T> &x,
516 const Array<OneD, const int> &y)
517{
518 ASSERTL1(n <= w.size() + w.GetOffset(), "Array out of bounds");
519 ASSERTL1(n <= x.size() + x.GetOffset(), "Array out of bounds");
520 ASSERTL1(n <= y.size() + y.GetOffset(), "Array out of bounds");
521
522 return Dot2(n, &w[0], &x[0], &y[0]);
523}
524
525/// \brief dot product
526template <class T>
527T Ddot(int n, const Array<OneD, const T> &w, const int incw,
528 const Array<OneD, const T> &x, const int incx,
529 const Array<OneD, const int> &y, const int incy)
530{
531 ASSERTL1(n * incw <= w.size() + w.GetOffset(), "Array out of bounds");
532 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
533 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
534
535 return Dot2(n, &w[0], incw, &x[0], incx, &y[0], incy);
536}
537
538/********** Memory routines ***********************/
539
540// \brief copy one vector to another
541template <class T>
542void Vcopy(int n, const Array<OneD, const T> &x, int incx, Array<OneD, T> &y,
543 int const incy)
544{
545 ASSERTL1(static_cast<unsigned int>(std::abs(n * incx)) <=
546 x.size() + x.GetOffset(),
547 "Array out of bounds");
548 ASSERTL1(static_cast<unsigned int>(std::abs(n * incy)) <=
549 y.size() + y.GetOffset(),
550 "Array out of bounds");
551
552 Vcopy(n, &x[0], incx, &y[0], incy);
553}
554
555// \brief reverse the ordering of vector to another
556template <class T>
557void Reverse(int n, const Array<OneD, const T> &x, int incx, Array<OneD, T> &y,
558 int const incy)
559{
560 ASSERTL1(static_cast<unsigned int>(std::abs(n * incx)) <=
561 x.size() + x.GetOffset(),
562 "Array out of bounds");
563 ASSERTL1(static_cast<unsigned int>(std::abs(n * incy)) <=
564 y.size() + y.GetOffset(),
565 "Array out of bounds");
566
567 Reverse(n, &x[0], incx, &y[0], incy);
568}
569
570} // namespace Vmath
571#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:250
Nektar::LibUtilities::beta
@ beta
Gauss Radau pinned at x=-1,.
Definition PointsType.h:59
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:527
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
Generated on Sat Jul 4 2026 01:02:21 for Nektar++ by doxygen 1.9.8