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 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Collection of templated functions for vector mathematics
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #ifndef NEKTAR_LIB_LIBUTILITIES_BASSICUTILS_VECTORMATH_HPP
37 #define NEKTAR_LIB_LIBUTILITIES_BASSICUTILS_VECTORMATH_HPP
38 
41 
42 namespace Vmath
43 {
44 
45  /***************** Math routines ***************/
46 
47  /// \brief Fill a vector with a constant value
48  template<class T> LIB_UTILITIES_EXPORT void Fill( int n, const T alpha, T *x, const int incx );
49 
50 
51 
52  /// \brief Generates a number from ~Normal(0,1)
53  template<class T> LIB_UTILITIES_EXPORT T ran2 (long* idum);
54 
55 
56  /// \brief Fills a vector with white noise.
57  template<class T> LIB_UTILITIES_EXPORT void FillWhiteNoise( int n, const T eps, T *x, const int incx, int seed = 0);
58 
59  /// \brief Multiply vector z = x*y
60  template<class T> LIB_UTILITIES_EXPORT void Vmul( int n, const T *x, const int incx, const T *y,
61  const int incy, T*z, const int incz);
62 
63  /// \brief Scalar multiply y = alpha*y
64 
65  template<class T> LIB_UTILITIES_EXPORT void Smul( int n, const T alpha, const T *x, const int incx,
66  T *y, const int incy);
67 
68  /// \brief Multiply vector z = x/y
69  template<class T> LIB_UTILITIES_EXPORT void Vdiv( int n, const T *x, const int incx, const T *y,
70  const int incy, T*z, const int incz);
71 
72  /// \brief Scalar multiply y = alpha/y
73  template<class T> LIB_UTILITIES_EXPORT void Sdiv( int n, const T alpha, const T *x,
74  const int incx, T *y, const int incy);
75 
76  /// \brief Add vector z = x+y
77  template<class T> LIB_UTILITIES_EXPORT void Vadd( int n, const T *x, const int incx, const T *y,
78  const int incy, T *z, const int incz);
79 
80  /// \brief Add vector y = alpha + x
81  template<class T> LIB_UTILITIES_EXPORT void Sadd( int n, const T alpha, const T *x,
82  const int incx, T *y, const int incy);
83 
84  /// \brief Subtract vector z = x-y
85  template<class T> LIB_UTILITIES_EXPORT void Vsub( int n, const T *x, const int incx, const T *y,
86  const int incy, T *z, const int incz);
87 
88  /// \brief Zero vector
89  template<class T> LIB_UTILITIES_EXPORT void Zero(int n, T *x, const int incx);
90 
91  /// \brief Negate x = -x
92  template<class T> LIB_UTILITIES_EXPORT void Neg( int n, T *x, const int incx);
93 
94 
95  template<class T> void Vlog(int n, const T *x, const int incx,
96  T *y, const int incy)
97  {
98  while (n--)
99  {
100  *y = log( *x );
101  x += incx;
102  y += incy;
103  }
104  }
105 
106 
107  template<class T> void Vexp(int n, const T *x, const int incx,
108  T *y, const int incy)
109  {
110  while (n--)
111  {
112  *y = exp( *x );
113  x += incx;
114  y += incy;
115  }
116  }
117 
118  template<class T> void Vpow(int n, const T *x, const int incx,
119  const T f, T *y, const int incy)
120  {
121  while (n--)
122  {
123  *y = pow( *x, f );
124  x += incx;
125  y += incy;
126  }
127  }
128 
129 
130  /// \brief sqrt y = sqrt(x)
131  template<class T> LIB_UTILITIES_EXPORT void Vsqrt(int n, const T *x, const int incx,
132  T *y, const int incy);
133 
134  /// \brief vabs: y = |x|
135  template<class T> LIB_UTILITIES_EXPORT void Vabs(int n, const T *x, const int incx,
136  T *y, const int incy);
137 
138  /********** Triad routines ***********************/
139 
140  /// \brief vvtvp (vector times vector plus vector): z = w*x + y
141  template<class T> LIB_UTILITIES_EXPORT void Vvtvp(int n,
142  const T *w, const int incw,
143  const T *x, const int incx,
144  const T *y, const int incy,
145  T *z, const int incz);
146 
147  /// \brief vvtvm (vector times vector plus vector): z = w*x - y
148  template<class T> LIB_UTILITIES_EXPORT void Vvtvm(int n, const T *w, const int incw, const T *x,
149  const int incx, const T *y, const int incy,
150  T *z, const int incz);
151 
152  /// \brief svtvp (scalar times vector plus vector): z = alpha*x + y
153  template<class T> LIB_UTILITIES_EXPORT void Svtvp(int n, const T alpha, const T *x,
154  const int incx, const T *y, const int incy,
155  T *z, const int incz);
156 
157  /// \brief svtvp (scalar times vector plus vector): z = alpha*x - y
158  template<class T> LIB_UTILITIES_EXPORT void Svtvm(int n, const T alpha, const T *x,
159  const int incx, const T *y, const int incy,
160  T *z, const int incz);
161 
162  /// \brief vvtvvtp (vector times vector plus vector times vector):
163  // z = v*w + x*y
164  template<class T> LIB_UTILITIES_EXPORT void Vvtvvtp (int n,
165  const T* v, int incv,
166  const T* w, int incw,
167  const T* x, int incx,
168  const T* y, int incy,
169  T* z, int incz);
170  /// \brief vvtvvtm (vector times vector minus vector times vector):
171  // z = v*w - x*y
172  template<class T> LIB_UTILITIES_EXPORT void Vvtvvtm (int n,
173  const T* v, int incv,
174  const T* w, int incw,
175  const T* x, int incx,
176  const T* y, int incy,
177  T* z, int incz);
178  /// \brief Svtsvtp (scalar times vector plus scalar times vector):
179  // z = alpha*x + beta*y
180  template<class T> LIB_UTILITIES_EXPORT void Svtsvtp (int n,
181  const T alpha,
182  const T* x, int incx,
183  const T beta,
184  const T* y, int incy,
185  T* z, int incz);
186 
187  /// \brief Vstvpp (scalar times vector plus vector plus vector):
188  // z = v*w + x*y
189  template<class T> LIB_UTILITIES_EXPORT void Vstvpp(int n,
190  const T alpha,
191  const T* v, int incv,
192  const T* w, int incw,
193  const T* x, int incx,
194  T* z, int incz);
195 
196  /************ Misc routine from Veclib (and extras) ************/
197 
198  /// \brief Gather vector z[i] = x[y[i]]
199  template<class T> LIB_UTILITIES_EXPORT void Gathr(int n, const T *x, const int *y,
200  T *z);
201 
202  /// \brief Gather vector z[i] = sign[i]*x[y[i]]
203  template<class T> LIB_UTILITIES_EXPORT void Gathr(int n, const T *sign, const T *x, const int *y,
204  T *z);
205 
206  /// \brief Scatter vector z[y[i]] = x[i]
207  template<class T> LIB_UTILITIES_EXPORT void Scatr(int n, const T *x, const int *y,
208  T *z);
209 
210  /// \brief Scatter vector z[y[i]] = sign[i]*x[i]
211  template<class T> LIB_UTILITIES_EXPORT void Scatr(int n, const T *sign, const T *x, const int *y,
212  T *z);
213 
214  /// \brief Assemble z[y[i]] += x[i]; z should be zero'd first
215  template<class T> LIB_UTILITIES_EXPORT void Assmb(int n, const T *x, const int *y,
216  T *z);
217 
218  /// \brief Assemble z[y[i]] += sign[i]*x[i]; z should be zero'd first
219  template<class T> LIB_UTILITIES_EXPORT void Assmb(int n, const T *sign, const T *x, const int *y,
220  T *z);
221 
222 
223  /************* Reduction routines *****************/
224 
225  /// \brief Subtract return sum(x)
226  template<class T> LIB_UTILITIES_EXPORT T Vsum( int n, const T *x, const int incx);
227 
228 
229  /// \brief Return the index of the maximum element in x
230  template<class T> LIB_UTILITIES_EXPORT int Imax( int n, const T *x, const int incx);
231 
232  /// \brief Return the maximum element in x -- called vmax to avoid
233  /// conflict with max
234  template<class T> LIB_UTILITIES_EXPORT T Vmax( int n, const T *x, const int incx);
235 
236  /// \brief Return the index of the maximum absolute element in x
237  template<class T> LIB_UTILITIES_EXPORT int Iamax( int n, const T *x, const int incx);
238 
239  /// \brief Return the maximum absolute element in x
240  /// called vamax to avoid conflict with max
241  template<class T> LIB_UTILITIES_EXPORT T Vamax( int n, const T *x, const int incx);
242 
243 
244  /// \brief Return the index of the minimum element in x
245  template<class T> LIB_UTILITIES_EXPORT int Imin( int n, const T *x, const int incx);
246 
247 
248  /// \brief Return the minimum element in x - called vmin to avoid
249  /// conflict with min
250  template<class T> LIB_UTILITIES_EXPORT T Vmin( int n, const T *x, const int incx);
251 
252  /// \brief Return number of NaN elements of x
253  template<class T> LIB_UTILITIES_EXPORT int Nnan(int n, const T *x, const int incx);
254 
255 
256  /// \brief vvtvp (vector times vector times vector): z = w*x*y
257  template<class T> LIB_UTILITIES_EXPORT T Dot( int n,
258  const T *w,
259  const T *x);
260 
261  /// \brief vvtvp (vector times vector times vector): z = w*x*y
262  template<class T> LIB_UTILITIES_EXPORT T Dot( int n,
263  const T *w, const int incw,
264  const T *x, const int incx);
265 
266  /// \brief vvtvp (vector times vector times vector): z = w*x*y
267  template<class T> LIB_UTILITIES_EXPORT T Dot2( int n,
268  const T *w,
269  const T *x,
270  const int *y);
271 
272  /// \brief vvtvp (vector times vector times vector): z = w*x*y
273  template<class T> LIB_UTILITIES_EXPORT T Dot2( int n,
274  const T *w, const int incw,
275  const T *x, const int incx,
276  const int *y, const int incy);
277 
278  /********** Memory routines ***********************/
279 
280  // \brief copy one vector to another
281  template<class T>
282  LIB_UTILITIES_EXPORT void Vcopy(int n, const T *x, const int incx,
283  T *y, const int incy);
284 
285  // \brief reverse the ordering of vector to another
286  template<class T> LIB_UTILITIES_EXPORT void Reverse( int n, const T *x, const int incx, T *y, const int incy);
287 
288 
289 }
290 #endif //VECTORMATH_HPP
291 
Definition: Vmath.cpp:40
void Vpow(int n, const T *x, const int incx, const T f, T *y, const int incy)
Definition: Vmath.hpp:118
void Gathr(int n, const T *x, const int *y, T *z)
Gather vector z[i] = x[y[i]].
Definition: Vmath.cpp:630
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:22
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:394
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:756
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:848
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
Definition: Vmath.cpp:824
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:471
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:428
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/y.
Definition: Vmath.cpp:257
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:227
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
Definition: Vmath.cpp:410
int Imax(int n, const T *x, const int incx)
Return the index of the maximum element in x.
Definition: Vmath.cpp:732
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1062
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
T ran2(long *idum)
Generates a number from ~Normal(0,1)
Definition: Vmath.cpp:73
#define LIB_UTILITIES_EXPORT
int Nnan(int n, const T *x, const int incx)
Return number of NaN elements of x.
Definition: Vmath.cpp:869
void Vexp(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:107
void Svtvm(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x - y
Definition: Vmath.cpp:504
void Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
Definition: Vmath.cpp:659
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
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:686
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:550
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:802
int Iamax(int n, const T *x, const int incx)
Return the index of the maximum absolute element in x.
Definition: Vmath.cpp:777
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.cpp:301
T Dot(int n, const T *w, const T *x)
vvtvp (vector times vector times vector): z = w*x*y
Definition: Vmath.cpp:891
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:329
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:603
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:523
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 plus vector): z = w*x - y
Definition: Vmath.cpp:451
T Dot2(int n, const T *w, const T *x, const int *y)
vvtvp (vector times vector times vector): z = w*x*y
Definition: Vmath.cpp:931
void FillWhiteNoise(int n, const T eps, T *x, const int incx, int outseed)
Fills a vector with white noise.
Definition: Vmath.cpp:138
void Vlog(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:95
void Svtsvtp(int n, const T alpha, const T *x, int incx, const T beta, const T *y, int incy, T *z, int incz)
vvtvvtp (scalar times vector plus scalar times vector):
Definition: Vmath.cpp:577
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:714
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
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:285
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:169