Nektar++
Vmath.hpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: Vmath.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: Collection of templated functions for vector mathematics
32//
33// Note: For those unfamiliar with the vector routines notation, it is
34// reverse polish notation (RPN). For example:
35//
36// In the function "Vvtvp()", it is "Vvt" means vector vector times,
37// which in infix notation is "v * v". So "Vvtvp" is:
38//
39// RPN: vector vector times vector plus
40// Infix: ( vector * vector ) + vector
41//
42///////////////////////////////////////////////////////////////////////////////
43
44#ifndef NEKTAR_LIB_LIBUTILITIES_BASSICUTILS_VECTORMATH_HPP
45#define NEKTAR_LIB_LIBUTILITIES_BASSICUTILS_VECTORMATH_HPP
46
47#include <LibUtilities/LibUtilitiesDeclspec.h>
48#include <LibUtilities/LinearAlgebra/Blas.hpp>
49
50namespace Vmath
51{
52/***************** Math routines ***************/
53
54/// \brief Fill a vector with a constant value
55template <class T>
56LIB_UTILITIES_EXPORT void Fill(int n, const T alpha, T *x, const int incx);
57
58/// \brief Generates a number from ~Normal(0,1)
59template <class T> LIB_UTILITIES_EXPORT T ran2(long *idum);
60
61/// \brief Fills a vector with white noise.
62template <class T>
63LIB_UTILITIES_EXPORT void FillWhiteNoise(int n, const T eps, T *x,
64 const int incx, int seed = 9999);
65
66/// \brief Multiply vector z = x*y
67template <class T>
68LIB_UTILITIES_EXPORT void Vmul(int n, const T *x, const int incx, const T *y,
69 const int incy, T *z, const int incz);
70
71/// \brief Scalar multiply y = alpha*x
72
73template <class T>
74LIB_UTILITIES_EXPORT void Smul(int n, const T alpha, const T *x, const int incx,
75 T *y, const int incy);
76
77/// \brief Multiply vector z = x/y
78template <class T>
79LIB_UTILITIES_EXPORT void Vdiv(int n, const T *x, const int incx, const T *y,
80 const int incy, T *z, const int incz);
81
82/// \brief Scalar multiply y = alpha/x
83template <class T>
84LIB_UTILITIES_EXPORT void Sdiv(int n, const T alpha, const T *x, const int incx,
85 T *y, const int incy);
86
87/// \brief Add vector z = x+y
88template <class T>
89LIB_UTILITIES_EXPORT void Vadd(int n, const T *x, const int incx, const T *y,
90 const int incy, T *z, const int incz);
91
92/// \brief Add vector y = alpha + x
93template <class T>
94LIB_UTILITIES_EXPORT void Sadd(int n, const T alpha, const T *x, const int incx,
95 T *y, const int incy);
96
97/// \brief Subtract vector z = x-y
98template <class T>
99LIB_UTILITIES_EXPORT void Vsub(int n, const T *x, const int incx, const T *y,
100 const int incy, T *z, const int incz);
101
102/// \brief Substract vector y = alpha - x
103template <class T>
104LIB_UTILITIES_EXPORT void Ssub(int n, const T alpha, const T *x, const int incx,
105 T *y, const int incy);
106
107/// \brief Zero vector
108template <class T> LIB_UTILITIES_EXPORT void Zero(int n, T *x, const int incx);
109
110/// \brief Negate x = -x
111template <class T> LIB_UTILITIES_EXPORT void Neg(int n, T *x, const int incx);
112
113/// \brief log y = log(x)
114template <class T>
115void Vlog(int n, const T *x, const int incx, T *y, const int incy)
116{
117 while (n--)
118 {
119 *y = log(*x);
120 x += incx;
121 y += incy;
122 }
123}
124
125/// \brief exp y = exp(x)
126template <class T>
127void Vexp(int n, const T *x, const int incx, T *y, const int incy)
128{
129 while (n--)
130 {
131 *y = exp(*x);
132 x += incx;
133 y += incy;
134 }
135}
136
137/// \brief pow y = pow(x, f)
138template <class T>
139void Vpow(int n, const T *x, const int incx, const T f, T *y, const int incy)
140{
141 while (n--)
142 {
143 *y = pow(*x, f);
144 x += incx;
145 y += incy;
146 }
147}
148
149/// \brief sqrt y = sqrt(x)
150template <class T>
151LIB_UTILITIES_EXPORT void Vsqrt(int n, const T *x, const int incx, T *y,
152 const int incy);
153
154/// \brief vabs: y = |x|
155template <class T>
156LIB_UTILITIES_EXPORT void Vabs(int n, const T *x, const int incx, T *y,
157 const int incy);
158
159/********** Triad routines ***********************/
160
161/// \brief vvtvp (vector times vector plus vector): z = w*x + y
162template <class T>
163LIB_UTILITIES_EXPORT void Vvtvp(int n, const T *w, const int incw, const T *x,
164 const int incx, const T *y, const int incy,
165 T *z, const int incz);
166
167/// \brief vvtvm (vector times vector minus vector): z = w*x - y
168template <class T>
169LIB_UTILITIES_EXPORT void Vvtvm(int n, const T *w, const int incw, const T *x,
170 const int incx, const T *y, const int incy,
171 T *z, const int incz);
172
173/// \brief Svtvp (scalar times vector plus vector): z = alpha*x + y
174template <class T>
175LIB_UTILITIES_EXPORT void Svtvp(int n, const T alpha, const T *x,
176 const int incx, const T *y, const int incy,
177 T *z, const int incz);
178
179/// \brief Svtvm (scalar times vector minus vector): z = alpha*x - y
180template <class T>
181LIB_UTILITIES_EXPORT void Svtvm(int n, const T alpha, const T *x,
182 const int incx, const T *y, const int incy,
183 T *z, const int incz);
184
185/// \brief vvtvvtp (vector times vector plus vector times vector):
186// z = v*w + x*y
187template <class T>
188LIB_UTILITIES_EXPORT void Vvtvvtp(int n, const T *v, int incv, const T *w,
189 int incw, const T *x, int incx, const T *y,
190 int incy, T *z, int incz);
191
192/// \brief vvtvvtm (vector times vector minus vector times vector):
193// z = v*w - x*y
194template <class T>
195LIB_UTILITIES_EXPORT void Vvtvvtm(int n, const T *v, int incv, const T *w,
196 int incw, const T *x, int incx, const T *y,
197 int incy, T *z, int incz);
198
199/// \brief Svtsvtp (scalar times vector plus scalar times vector):
200// z = alpha*x + beta*y
201template <class T>
202LIB_UTILITIES_EXPORT void Svtsvtp(int n, const T alpha, const T *x, int incx,
203 const T beta, const T *y, int incy, T *z,
204 int incz);
205
206/// \brief Vstvpp (scalar times vector plus vector plus vector):
207// z = alpha*w + x + y
208template <class T>
209LIB_UTILITIES_EXPORT void Vstvpp(int n, const T alpha, const T *v, int incv,
210 const T *w, int incw, const T *x, int incx,
211 T *z, int incz);
212
213/************ Misc routine from Veclib (and extras) ************/
214
215/// \brief Gather vector z[i] = x[y[i]]
216template <class T, class I,
217 typename = typename std::enable_if<std::is_floating_point<T>::value &&
218 std::is_integral<I>::value>::type>
219void Gathr(I n, const T *x, const I *y, T *z)
220{
221 while (n--)
222 {
223 *z++ = *(x + *y++);
224 }
225 return;
226}
227
228/// \brief Gather vector z[i] = sign[i]*x[y[i]]
229template <class T>
230LIB_UTILITIES_EXPORT void Gathr(int n, const T *sign, const T *x, const int *y,
231 T *z);
232
233/// \brief Scatter vector z[y[i]] = x[i]
234template <class T>
235LIB_UTILITIES_EXPORT void Scatr(int n, const T *x, const int *y, T *z);
236
237/// \brief Scatter vector z[y[i]] = sign[i]*x[i]
238template <class T>
239LIB_UTILITIES_EXPORT void Scatr(int n, const T *sign, const T *x, const int *y,
240 T *z);
241
242/// \brief Assemble z[y[i]] += x[i]; z should be zero'd first
243template <class T>
244LIB_UTILITIES_EXPORT void Assmb(int n, const T *x, const int *y, T *z);
245
246/// \brief Assemble z[y[i]] += sign[i]*x[i]; z should be zero'd first
247template <class T>
248LIB_UTILITIES_EXPORT void Assmb(int n, const T *sign, const T *x, const int *y,
249 T *z);
250
251/************* Reduction routines *****************/
252
253/// \brief Subtract return sum(x)
254template <class T>
255LIB_UTILITIES_EXPORT T Vsum(int n, const T *x, const int incx);
256
257/// \brief Return the index of the maximum element in x
258template <class T>
259LIB_UTILITIES_EXPORT int Imax(int n, const T *x, const int incx);
260
261/// \brief Return the maximum element in x -- called vmax to avoid
262/// conflict with max
263template <class T>
264LIB_UTILITIES_EXPORT T Vmax(int n, const T *x, const int incx);
265
266/// \brief Return the index of the maximum absolute element in x
267template <class T>
268LIB_UTILITIES_EXPORT int Iamax(int n, const T *x, const int incx);
269
270/// \brief Return the maximum absolute element in x
271/// called vamax to avoid conflict with max
272template <class T>
273LIB_UTILITIES_EXPORT T Vamax(int n, const T *x, const int incx);
274
275/// \brief Return the index of the minimum element in x
276template <class T>
277LIB_UTILITIES_EXPORT int Imin(int n, const T *x, const int incx);
278
279/// \brief Return the minimum element in x - called vmin to avoid
280/// conflict with min
281template <class T>
282LIB_UTILITIES_EXPORT T Vmin(int n, const T *x, const int incx);
283
284/// \brief Return number of NaN elements of x
285template <class T>
286LIB_UTILITIES_EXPORT int Nnan(int n, const T *x, const int incx);
287
288/// \brief dot product
289template <class T> LIB_UTILITIES_EXPORT T Dot(int n, const T *w, const T *x);
290
291/// \brief dot product
292template <class T>
293LIB_UTILITIES_EXPORT T Dot(int n, const T *w, const int incw, const T *x,
294 const int incx);
295
296/// \brief dot product
297template <class T>
298LIB_UTILITIES_EXPORT T Dot2(int n, const T *w, const T *x, const int *y);
299
300/// \brief dot product
301template <class T>
302LIB_UTILITIES_EXPORT T Dot2(int n, const T *w, const int incw, const T *x,
303 const int incx, const int *y, const int incy);
304
305/********** Memory routines ***********************/
306
307// \brief copy one vector to another
308template <class T>
309LIB_UTILITIES_EXPORT void Vcopy(int n, const T *x, const int incx, T *y,
310 const int incy);
311
312// \brief reverse the ordering of vector to another
313template <class T>
314LIB_UTILITIES_EXPORT void Reverse(int n, const T *x, const int incx, T *y,
315 const int incy);
316
317} // namespace Vmath
318#endif // VECTORMATH_HPP
LIB_UTILITIES_EXPORT
#define LIB_UTILITIES_EXPORT
Definition: LibUtilitiesDeclspec.h:45
sign
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:49
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)
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::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::ran2
T ran2(long *idum)
Generates a number from ~Normal(0,1)
Definition: Vmath.cpp:73
Vmath::Vstvpp
void Vstvpp(int n, const T alpha, const T *v, int incv, const T *w, int incw, const T *x, int incx, T *z, int incz)
Vstvpp (scalar times vector plus vector plus vector):
Definition: Vmath.cpp:772
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::log
scalarT< T > log(scalarT< T > in)
Definition: scalar.hpp:303
Generated on Wed Jan 31 2024 08:26:46 for Nektar++ by doxygen 1.9.4