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, 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
87template <class T>
88void 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*x
100
101template <class T>
102void 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
114template <class T>
115void 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/x
130template <class T>
131void 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
143template <class T>
144void 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
167template <class T>
168void 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
179template <class T>
180void 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
192template <class T>
193void 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
204template <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
212template <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/// \brief log y = log(x)
220template <class T>
221void Vlog(int n, const Array<OneD, const T> &x, const int incx,
222 Array<OneD, T> &y, const int incy)
223{
224 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
225 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
226
227 Vlog(n, &x[0], incx, &y[0], incy);
228}
229
230/// \brief exp y = exp(x)
231template <class T>
232void Vexp(int n, const Array<OneD, const T> &x, const int incx,
233 Array<OneD, T> &y, const int incy)
234{
235 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
236 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
237
238 Vexp(n, &x[0], incx, &y[0], incy);
239}
240
241/// \brief pow y = pow(x, f)
242template <class T>
243void Vpow(int n, const Array<OneD, const T> &x, const int incx, const T f,
244 Array<OneD, T> &y, const int incy)
245{
246 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
247 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
248
249 Vpow(n, &x[0], incx, f, &y[0], incy);
250}
251
252/// \brief sqrt y = sqrt(x)
253template <class T>
254void Vsqrt(int n, const Array<OneD, const T> &x, const int incx,
255 Array<OneD, T> &y, const int incy)
256{
257 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
258 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
259
260 Vsqrt(n, &x[0], incx, &y[0], incy);
261}
262
263/// \brief vabs: y = |x|
264template <class T>
265void Vabs(int n, const Array<OneD, const T> &x, const int incx,
266 Array<OneD, T> &y, const int incy)
267{
268 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
269 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
270
271 Vabs(n, &x[0], incx, &y[0], incy);
272}
273
274/********** Triad routines ***********************/
275
276/// \brief vvtvp (vector times vector plus vector): z = w*x + y
277template <class T>
278void Vvtvp(int n, const Array<OneD, const T> &w, const int incw,
279 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() + w.GetOffset(), "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#ifdef NEKTAR_ENABLE_SIMD_VMATH
289 boost::ignore_unused(incw, incx, incy, incz);
290 ASSERTL1(incw == 1, "Simd vmath requires inc = 1");
291 ASSERTL1(incx == 1, "Simd vmath requires inc = 1");
292 ASSERTL1(incy == 1, "Simd vmath requires inc = 1");
293 ASSERTL1(incz == 1, "Simd vmath requires inc = 1");
294 SIMD::Vvtvp(n, &w[0], &x[0], &y[0], &z[0]);
295#else
296 Vvtvp(n, &w[0], incw, &x[0], incx, &y[0], incy, &z[0], incz);
297#endif
298}
299
300/// \brief Vvtvp (vector times vector plus vector): z = w*x + y
301template <class T>
302void Vvtvp(int n, const Array<TwoD, NekDouble>::const_reference &w,
303 const int incw, const Array<OneD, const T> &x, const int incx,
304 const Array<OneD, const T> &y, const int incy, Array<OneD, T> &z,
305 const int incz)
306{
307 ASSERTL1(n * incw <= w.size(), "Array out of bounds");
308 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
309 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
310 ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
311
312 Vvtvp(n, w.origin(), incw, &x[0], incx, &y[0], incy, &z[0], incz);
313}
314
315/// \brief vvtvm (vector times vector minus vector): z = w*x - y
316template <class T>
317void Vvtvm(int n, const Array<OneD, const T> &w, const int incw,
318 const Array<OneD, const T> &x, const int incx,
319 const Array<OneD, const T> &y, const int incy, Array<OneD, T> &z,
320 const int incz)
321{
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(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::Vvtvm(n, &w[0], &x[0], &y[0], &z[0]);
334#else
335 Vvtvm(n, &w[0], incw, &x[0], incx, &y[0], incy, &z[0], incz);
336#endif
337}
338
339/// \brief svtvp (scalar times vector plus vector): z = alpha*x + y
340template <class T>
341void Svtvp(int n, const T alpha, const Array<OneD, const T> &x, const int incx,
342 const Array<OneD, const T> &y, const int incy, Array<OneD, T> &z,
343 const int incz)
344{
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 Svtvp(n, alpha, &x[0], incx, &y[0], incy, &z[0], incz);
350}
351
352/// \brief svtvm (scalar times vector minus vector): z = alpha*x - y
353template <class T>
354void Svtvm(int n, const T alpha, const Array<OneD, const T> &x, const int incx,
355 const Array<OneD, const T> &y, const int incy, Array<OneD, T> &z,
356 const int incz)
357{
358 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
359 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
360 ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
361
362 Svtvm(n, alpha, &x[0], incx, &y[0], incy, &z[0], incz);
363}
364
365/// \brief vvtvvtp (vector times vector plus vector times vector): z = v*w + x*y
366template <class T>
367void Vvtvvtp(int n, const Array<OneD, const T> &v, int incv,
368 const Array<OneD, const T> &w, int incw,
369 const Array<OneD, const T> &x, int incx,
370 const Array<OneD, const T> &y, int incy, Array<OneD, T> &z,
371 int incz)
372{
373 ASSERTL1(n * incv <= v.size() + v.GetOffset(), "Array out of bounds");
374 ASSERTL1(n * incw <= w.size() + w.GetOffset(), "Array out of bounds");
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#ifdef NEKTAR_ENABLE_SIMD_VMATH
380 boost::ignore_unused(incv, incw, incx, incy, incz);
381 ASSERTL1(incw == 1, "Simd vmath requires inc = 1");
382 ASSERTL1(incx == 1, "Simd vmath requires inc = 1");
383 ASSERTL1(incy == 1, "Simd vmath requires inc = 1");
384 ASSERTL1(incz == 1, "Simd vmath requires inc = 1");
385 SIMD::Vvtvvtp(n, &v[0], &w[0], &x[0], &y[0], &z[0]);
386#else
387 Vvtvvtp(n, &v[0], incv, &w[0], incw, &x[0], incx, &y[0], incy, &z[0], incz);
388#endif
389}
390
391/// \brief vvtvvtm (vector times vector minus vector times vector): z = v*w -
392/// x*y
393template <class T>
394void Vvtvvtm(int n, const Array<OneD, const T> &v, int incv,
395 const Array<OneD, const T> &w, int incw,
396 const Array<OneD, const T> &x, int incx,
397 const Array<OneD, const T> &y, int incy, Array<OneD, T> &z,
398 int incz)
399{
400 ASSERTL1(n * incv <= v.size() + v.GetOffset(), "Array out of bounds");
401 ASSERTL1(n * incw <= w.size() + w.GetOffset(), "Array out of bounds");
402 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
403 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
404 ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
405
406#ifdef NEKTAR_ENABLE_SIMD_VMATH
407 boost::ignore_unused(incv, incw, incx, incy, incz);
408 ASSERTL1(incw == 1, "Simd vmath requires inc = 1");
409 ASSERTL1(incx == 1, "Simd vmath requires inc = 1");
410 ASSERTL1(incy == 1, "Simd vmath requires inc = 1");
411 ASSERTL1(incz == 1, "Simd vmath requires inc = 1");
412 SIMD::Vvtvvtm(n, &v[0], &w[0], &x[0], &y[0], &z[0]);
413#else
414 Vvtvvtm(n, &v[0], incv, &w[0], incw, &x[0], incx, &y[0], incy, &z[0], incz);
415#endif
416}
417
418/// \brief svtsvtp (scalar times vector plus scalar times vector): z = alpha*x +
419/// beta*y
420template <class T>
421void Svtsvtp(int n, const T alpha, const Array<OneD, const T> &x,
422 const int incx, const T beta, const Array<OneD, const T> &y,
423 const int incy, Array<OneD, T> &z, const int incz)
424{
425 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
426 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
427 ASSERTL1(n * incz <= z.size() + z.GetOffset(), "Array out of bounds");
428
429 Svtsvtp(n, alpha, &x[0], incx, beta, &y[0], incy, &z[0], incz);
430}
431
432/************ Misc routine from Veclib (and extras) ************/
433
434/// \brief Gather vector z[i] = x[y[i]]
435template <class T, class I,
436 typename = typename std::enable_if<std::is_floating_point<T>::value &&
437 std::is_integral<I>::value>::type>
438void Gathr(I n, const Array<OneD, const T> &x, const Array<OneD, I> &y,
439 Array<OneD, T> &z)
440{
441 ASSERTL1(n <= y.size() + y.GetOffset(), "Array out of bounds");
442 ASSERTL1(n <= z.size() + z.GetOffset(), "Array out of bounds");
443
444#ifdef NEKTAR_ENABLE_SIMD_VMATH
445 SIMD::Gathr(n, &x[0], &y[0], &z[0]);
446#else
447 Gathr(n, &x[0], &y[0], &z[0]);
448#endif
449}
450
451/// \brief Scatter vector z[y[i]] = x[i]
452template <class T>
453void Scatr(int n, const Array<OneD, const T> &x,
454 const Array<OneD, const int> &y, Array<OneD, T> &z)
455{
456 ASSERTL1(n <= x.size() + x.GetOffset(), "Array out of bounds");
457 ASSERTL1(n <= y.size() + y.GetOffset(), "Array out of bounds");
458
459 Scatr(n, &x[0], &y[0], &z[0]);
460}
461
462/// \brief Assemble z[y[i]] += x[i]; z should be zero'd first
463template <class T>
464void Assmb(int n, const Array<OneD, T> &x, const Array<OneD, int> &y,
465 Array<OneD, T> &z)
466{
467 ASSERTL1(n <= x.size() + x.GetOffset(), "Array out of bounds");
468 ASSERTL1(n <= y.size() + y.GetOffset(), "Array out of bounds");
469
470 Assmb(n, &x[0], &y[0], &z[0]);
471}
472
473/************* Reduction routines *****************/
474
475/// \brief Subtract return sum(x)
476template <class T> T Vsum(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 Vsum(n, &x[0], incx);
481}
482
483/// \brief Return the index of the maximum element in x
484template <class T>
485int Imax(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 Imax(n, &x[0], incx);
490}
491
492/// \brief Return the maximum element in x -- called vmax to avoid
493/// conflict with max
494template <class T> T Vmax(int n, const Array<OneD, const T> &x, const int incx)
495{
496 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
497
498 return Vmax(n, &x[0], incx);
499}
500
501/// \brief Return the index of the maximum absolute element in x
502template <class T>
503int Iamax(int n, const Array<OneD, const T> &x, const int incx)
504{
505 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
506
507 return Iamax(n, &x[0], incx);
508}
509
510/// \brief Return the maximum absolute element in x
511/// called vamax to avoid conflict with max
512template <class T> T Vamax(int n, const Array<OneD, const T> &x, const int incx)
513{
514 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
515
516 return Vamax(n, &x[0], incx);
517}
518
519/// \brief Return the index of the minimum element in x
520template <class T>
521int Imin(int n, const Array<OneD, const T> &x, const int incx)
522{
523 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
524
525 return Imin(n, &x[0], incx);
526}
527
528/// \brief Return the minimum element in x - called vmin to avoid
529/// conflict with min
530template <class T> T Vmin(int n, const Array<OneD, const T> &x, const int incx)
531{
532 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
533
534 return Vmin(n, &x[0], incx);
535}
536
537/// \brief Return number of NaN elements of x
538template <class T>
539int Nnan(int n, const Array<OneD, const T> &x, const int incx)
540{
541 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
542
543 return Nnan(n, &x[0], incx);
544}
545
546/// \brief dot product
547template <class T>
548T Dot(int n, const Array<OneD, const T> &w, const Array<OneD, const T> &x)
549{
550 ASSERTL1(n <= w.size() + w.GetOffset(), "Array out of bounds");
551 ASSERTL1(n <= x.size() + x.GetOffset(), "Array out of bounds");
552
553 return Dot(n, &w[0], &x[0]);
554}
555
556/// \brief dot product
557template <class T>
558T Dot(int n, const Array<OneD, const T> &w, const int incw,
559 const Array<OneD, const T> &x, const int incx)
560{
561 ASSERTL1(n * incw <= w.size() + w.GetOffset(), "Array out of bounds");
562 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
563
564 return Dot(n, &w[0], incw, &x[0], incx);
565}
566
567/// \brief dot product
568template <class T>
569T Dot2(int n, const Array<OneD, const T> &w, const Array<OneD, const T> &x,
570 const Array<OneD, const int> &y)
571{
572 ASSERTL1(n <= w.size() + w.GetOffset(), "Array out of bounds");
573 ASSERTL1(n <= x.size() + x.GetOffset(), "Array out of bounds");
574 ASSERTL1(n <= y.size() + y.GetOffset(), "Array out of bounds");
575
576 return Dot2(n, &w[0], &x[0], &y[0]);
577}
578
579/// \brief dot product
580template <class T>
581T Ddot(int n, const Array<OneD, const T> &w, const int incw,
582 const Array<OneD, const T> &x, const int incx,
583 const Array<OneD, const int> &y, const int incy)
584{
585 ASSERTL1(n * incw <= w.size() + w.GetOffset(), "Array out of bounds");
586 ASSERTL1(n * incx <= x.size() + x.GetOffset(), "Array out of bounds");
587 ASSERTL1(n * incy <= y.size() + y.GetOffset(), "Array out of bounds");
588
589 return Dot2(n, &w[0], incw, &x[0], incx, &y[0], incy);
590}
591
592/********** Memory routines ***********************/
593
594// \brief copy one vector to another
595template <class T>
596void Vcopy(int n, const Array<OneD, const T> &x, int incx, Array<OneD, T> &y,
597 int const incy)
598{
599 ASSERTL1(static_cast<unsigned int>(std::abs(n * incx)) <=
600 x.size() + x.GetOffset(),
601 "Array out of bounds");
602 ASSERTL1(static_cast<unsigned int>(std::abs(n * incy)) <=
603 y.size() + y.GetOffset(),
604 "Array out of bounds");
605
606 Vcopy(n, &x[0], incx, &y[0], incy);
607}
608
609// \brief reverse the ordering of vector to another
610template <class T>
611void Reverse(int n, const Array<OneD, const T> &x, int incx, Array<OneD, T> &y,
612 int const incy)
613{
614 ASSERTL1(static_cast<unsigned int>(std::abs(n * incx)) <=
615 x.size() + x.GetOffset(),
616 "Array out of bounds");
617 ASSERTL1(static_cast<unsigned int>(std::abs(n * incy)) <=
618 y.size() + y.GetOffset(),
619 "Array out of bounds");
620
621 Reverse(n, &x[0], incx, &y[0], incy);
622}
623
624} // namespace Vmath
625#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::UnitTests::w
std::vector< double > w(NPUPPER)
Nektar::UnitTests::z
std::vector< double > z(NPUPPER)
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:261
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:49
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:309
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:416
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:474
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:358
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:155
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.cpp:529
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:453
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:746
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:207
Vmath::Vlog
void Vlog(int n, const T *x, const int incx, T *y, const int incy)
log y = log(x)
Definition: Vmath.hpp:115
Vmath::Vexp
void Vexp(int n, const T *x, const int incx, T *y, const int incy)
exp y = exp(x)
Definition: Vmath.hpp:127
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:617
Vmath::Vabs
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
Definition: Vmath.cpp:548
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:1142
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:513
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:1045
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:581
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:569
Vmath::Vsum
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:890
Vmath::Scatr
void Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
Definition: Vmath.cpp:817
Vmath::Dot
T Dot(int n, const T *w, const T *x)
dot (vector times vector): z = w*x
Definition: Vmath.cpp:1095
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:800
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:857
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.cpp:659
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:354
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:593
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:716
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:245
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.cpp:319
Vmath::Imax
int Imax(int n, const T *x, const int incx)
Return the index of the maximum element in x.
Definition: Vmath.cpp:913
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:280
Vmath::Imin
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
Definition: Vmath.cpp:1018
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:487
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:43
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:153
Vmath::Nnan
int Nnan(int n, const T *x, const int incx)
Return number of NaN elements of x.
Definition: Vmath.cpp:1071
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:994
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:379
Vmath::Reverse
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1222
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:940
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:687
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1191
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:414
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:139
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:966
tinysimd::abs
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:298
Generated on Wed Jan 31 2024 08:26:46 for Nektar++ by doxygen 1.9.4