Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 // 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: Wrappers around Vmath routines using Array<OneD,T> as arugments
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #ifndef NEKTAR_LIB_LIBUTILITIES_BASSICUTILS_VECTORMATHARRAY_HPP
37 #define NEKTAR_LIB_LIBUTILITIES_BASSICUTILS_VECTORMATHARRAY_HPP
38 
39 #include <LibUtilities/BasicUtils/SharedArray.hpp>
40 #include <LibUtilities/BasicUtils/Vmath.hpp>
41 
42  namespace Vmath
43  {
44  using namespace Nektar;
45 
46  /***************** Math routines ***************/
47  /// \brief Fill a vector with a constant value
48  template<class T> void Fill( int n, const T alpha, Array<OneD, T> &x, const int incx )
49  {
50 
51  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Out of bounds");
52 
53  Fill(n,alpha,&x[0],incx);
54  }
55 
56  template<class T> void FillWhiteNoise( int n, const T eps, Array<OneD, T> &x, const int incx, int outseed)
57  {
58  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Out of bounds");
59 
60  FillWhiteNoise(n,eps,&x[0],incx,outseed);
61  }
62 
63  /// \brief Multiply vector z = x*y
64  template<class T> void Vmul( int n, const Array<OneD, const T> &x, const int incx, const Array<OneD, const T> &y, const int incy, Array<OneD,T> &z, const int incz)
65  {
66  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
67  ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),"Array out of bounds");
68  ASSERTL1(n*incz <= z.num_elements()+z.GetOffset(),"Array out of bounds");
69 
70  Vmul(n,&x[0],incx,&y[0],incy,&z[0],incz);
71  }
72 
73  template<class T> void Vmul( int n, const Array<TwoD,NekDouble>::const_reference &x, const int incx, const Array<OneD,const T> &y, const int incy, Array<OneD,T> &z, const int incz)
74  {
75  ASSERTL1(n*incx <= x.num_elements(),"Array out of bounds");
76  ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),"Array out of bounds");
77  ASSERTL1(n*incz <= z.num_elements()+z.GetOffset(),"Array out of bounds");
78 
79  Vmul(n,x.origin(),incx,&y[0],incy,&z[0],incz);
80  }
81 
82  /// \brief Scalar multiply y = alpha*y
83 
84  template<class T> void Smul( int n, const T alpha, const Array<OneD,const T> &x, const int incx, Array<OneD,T> &y, const int incy)
85  {
86  ASSERTL1(static_cast<unsigned int>(n*incx) <= x.num_elements()+x.GetOffset(),"Array out of bounds");
87  ASSERTL1(static_cast<unsigned int>(n*incy) <= y.num_elements()+y.GetOffset(),"Array out of bounds");
88 
89  Smul(n,alpha, &x[0],incx,&y[0],incy);
90  }
91 
92  /// \brief Multiply vector z = x/y
93  template<class T> void Vdiv( int n, const Array<OneD,const T> &x, const int incx, const Array<OneD,const T> &y, const int incy, Array<OneD,T> &z, const int incz)
94  {
95  ASSERTL1(static_cast<unsigned int>(n*incx) <= x.num_elements()+x.GetOffset(),"Array out of bounds");
96  ASSERTL1(static_cast<unsigned int>(n*incy) <= y.num_elements()+y.GetOffset(),"Array out of bounds");
97  ASSERTL1(static_cast<unsigned int>(n*incz) <= z.num_elements()+z.GetOffset(),"Array out of bounds");
98 
99  Vdiv(n,&x[0],incx,&y[0],incy,&z[0],incz);
100 
101  }
102 
103  /// \brief Scalar multiply y = alpha/y
104  template<class T> void Sdiv( int n, const T alpha, const Array<OneD,const T> &x, const int incx, Array<OneD,T> &y, const int incy)
105  {
106  ASSERTL1(static_cast<unsigned int>(n*incx) <= x.num_elements()+x.GetOffset(),"Array out of bounds");
107  ASSERTL1(static_cast<unsigned int>(n*incy) <= y.num_elements()+y.GetOffset(),"Array out of bounds");
108 
109  Sdiv(n,alpha,&x[0],incx,&y[0],incy);
110  }
111 
112  /// \brief Add vector z = x+y
113  template<class T> void Vadd( int n, const Array<OneD,const T> &x, const int incx, const Array<OneD,const T> &y, const int incy, Array<OneD,T> &z, const int incz)
114  {
115 
116  ASSERTL1(static_cast<unsigned int>(n*incx) <= x.num_elements()+x.GetOffset(),"Array out of bounds");
117  ASSERTL1(static_cast<unsigned int>(n*incy) <= y.num_elements()+y.GetOffset(),"Array out of bounds");
118  ASSERTL1(static_cast<unsigned int>(n*incz) <= z.num_elements()+z.GetOffset(),"Array out of bounds");
119 
120  Vadd(n,&x[0],incx,&y[0],incy,&z[0],incz);
121  }
122 
123  /// \brief Add vector y = alpha + x
124  template<class T> void Sadd( int n, const T alpha, const Array<OneD,const T> &x,const int incx, Array<OneD,T> &y, const int incy)
125  {
126 
127  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
128  ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),"Array out of bounds");
129 
130  Sadd(n,alpha,&x[0],incx,&y[0],incy);
131  }
132 
133  /// \brief Subtract vector z = x-y
134  template<class T> void Vsub( int n, const Array<OneD,const T> &x, const int incx, const Array<OneD,const T> &y, const int incy, Array<OneD,T> &z, const int incz)
135  {
136  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
137  ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),"Array out of bounds");
138  ASSERTL1(n*incz <= z.num_elements()+z.GetOffset(),"Array out of bounds");
139 
140  Vsub(n,&x[0],incx,&y[0],incy,&z[0],incz);
141 
142  }
143 
144  /// \brief Zero vector
145  template<class T> void Zero(int n, Array<OneD,T> &x, const int incx)
146  {
147  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
148 
149  Zero(n,&x[0],incx);
150 
151  }
152 
153  /// \brief Negate x = -x
154  template<class T> void Neg( int n, Array<OneD,T> &x, const int incx)
155  {
156  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
157 
158  Neg(n,&x[0],incx);
159 
160  }
161 
162  template<class T> void Vlog(int n, const Array<OneD,const T> &x, const int incx, Array<OneD,T> &y, const int incy)
163  {
164  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
165  ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),"Array out of bounds");
166 
167  Vlog(n, &x[0], incx, &y[0], incy);
168  }
169 
170 
171  template<class T> void Vexp(int n, const Array<OneD,const T> &x, const int incx, Array<OneD,T> &y, const int incy)
172  {
173  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
174  ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),"Array out of bounds");
175 
176  Vexp(n, &x[0], incx, &y[0], incy);
177  }
178 
179  template<class T> void Vpow(int n, const Array<OneD,const T> &x, const int incx, const T f, Array<OneD,T> &y, const int incy)
180  {
181  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
182  ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),"Array out of bounds");
183 
184  Vpow(n, &x[0], incx, f, &y[0], incy);
185  }
186 
187  /// \brief sqrt y = sqrt(x)
188  template<class T> void Vsqrt(int n, const Array<OneD,const T> &x, const int incx, Array<OneD,T> &y, const int incy)
189  {
190  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
191  ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),"Array out of bounds");
192 
193  Vsqrt(n,&x[0],incx,&y[0],incy);
194  }
195 
196  /// \brief vabs: y = |x|
197  template<class T> void Vabs(int n, const Array<OneD,const T> &x, const int incx, Array<OneD,T> &y, const int incy)
198  {
199  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
200  ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),"Array out of bounds");
201 
202  Vabs(n,&x[0],incx,&y[0],incy);
203  }
204 
205  /********** Triad routines ***********************/
206 
207  /// \brief vvtvp (vector times vector plus vector): z = w*x + y
208  template<class T> void Vvtvp(int n, const Array<OneD, const T> &w, const int incw, const Array<OneD,const T> &x, const int incx, const Array<OneD, const T> &y, const int incy, Array<OneD,T> &z, const int incz)
209  {
210  ASSERTL1(n*incw <= w.num_elements()+w.GetOffset(),"Array out of bounds");
211  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
212  ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),"Array out of bounds");
213  ASSERTL1(n*incz <= z.num_elements()+z.GetOffset(),"Array out of bounds");
214 
215  Vvtvp(n,&w[0],incw,&x[0],incx,&y[0],incy,&z[0],incz);
216  }
217 
218  template<class T> void Vvtvp(int n, const Array<TwoD,NekDouble>::const_reference &w, const int incw, const Array<OneD, const T> &x, const int incx, const Array<OneD,const T> &y, const int incy, Array<OneD,T> &z, const int incz)
219  {
220  ASSERTL1(n*incw <= w.num_elements(),"Array out of bounds");
221  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
222  ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),"Array out of bounds");
223  ASSERTL1(n*incz <= z.num_elements()+z.GetOffset(),"Array out of bounds");
224 
225  Vvtvp(n,w.origin(),incw,&x[0],incx,&y[0],incy,&z[0],incz);
226  }
227 
228  /// \brief svtvp (scalar times vector plus vector): z = alpha*x + y
229  template<class T> void Svtvp(int n, const T alpha, const Array<OneD,const T> &x, const int incx, const Array<OneD, const T> &y, const int incy, Array<OneD,T> &z, const int incz)
230  {
231  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
232  ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),"Array out of bounds");
233  ASSERTL1(n*incz <= z.num_elements()+z.GetOffset(),"Array out of bounds");
234 
235  Svtvp(n,alpha,&x[0],incx,&y[0],incy,&z[0],incz);
236 
237  }
238 
239  /// \brief svtvp (scalar times vector plus vector): z = alpha*x + y
240  template<class T> void Svtvm(int n, const T alpha, const Array<OneD,const T> &x, const int incx, const Array<OneD, const T> &y, const int incy, Array<OneD,T> &z, const int incz)
241  {
242  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
243  ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),"Array out of bounds");
244  ASSERTL1(n*incz <= z.num_elements()+z.GetOffset(),"Array out of bounds");
245 
246  Svtvm(n,alpha,&x[0],incx,&y[0],incy,&z[0],incz);
247 
248  }
249 
250  /// \brief vvtvm (vector times vector minus vector): z = w*x - y
251  template<class T> void Vvtvm(int n, const Array<OneD,const T> &w, const int incw, const Array<OneD,const T> &x, const int incx, const Array<OneD,const T> &y, const int incy, Array<OneD,T> &z, const int incz)
252  {
253  ASSERTL1(n*incw <= w.num_elements()+w.GetOffset(),"Array out of bounds");
254  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
255  ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),"Array out of bounds");
256  ASSERTL1(n*incz <= z.num_elements()+z.GetOffset(),"Array out of bounds");
257 
258  Vvtvm(n,&w[0],incw,&x[0],incx,&y[0],incy,&z[0],incz);
259 
260  }
261 
262  /// \brief vvtvvtp (vector times vector plus vector times vector): z = v*w + y*z
263  template<class T> void Vvtvvtp (
264  int n,
265  const Array<OneD,const T> &v, int incv,
266  const Array<OneD,const T> &w, int incw,
267  const Array<OneD,const T> &x, int incx,
268  const Array<OneD,const T> &y, int incy,
269  Array<OneD, T> &z, int incz)
270  {
271  ASSERTL1(n*incv <= v.num_elements()+v.GetOffset(),"Array out of bounds");
272  ASSERTL1(n*incw <= w.num_elements()+w.GetOffset(),"Array out of bounds");
273  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
274  ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),"Array out of bounds");
275  ASSERTL1(n*incz <= z.num_elements()+z.GetOffset(),"Array out of bounds");
276 
277  Vvtvvtp(n,&v[0],incv,&w[0],incw,&x[0],incx,&y[0],incy,&z[0],incz);
278  }
279 
280  /// \brief svtsvtp (scalar times vector plus scalar times vector): z = alpha*x + beta*y
281  template<class T> void Svtsvtp(int n, const T alpha, const Array<OneD,const T> &x, const int incx, const T beta, const Array<OneD,const T> &y, const int incy, Array<OneD,T> &z, const int incz)
282  {
283  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
284  ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),"Array out of bounds");
285  ASSERTL1(n*incz <= z.num_elements()+z.GetOffset(),"Array out of bounds");
286 
287  Svtsvtp(n,alpha,&x[0],incx,beta,&y[0],incy,&z[0],incz);
288  }
289 
290 
291 
292 
293  /************ Misc routine from Veclib (and extras) ************/
294 
295  /// \brief Gather vector z[i] = x[y[i]]
296  template<class T> void Gathr(int n, const Array<OneD, const T> &x, const Array<OneD, const int> &y, Array<OneD,T> &z)
297  {
298 
299  ASSERTL1(n <= y.num_elements()+y.GetOffset(),"Array out of bounds");
300  ASSERTL1(n <= z.num_elements()+z.GetOffset(),"Array out of bounds");
301 
302  Gathr(n,&x[0],&y[0],&z[0]);
303 
304  }
305 
306  /// \brief Scatter vector z[y[i]] = x[i]
307  template<class T> void Scatr(int n, const Array<OneD,const T> &x, const Array<OneD,const int> &y, Array<OneD,T> &z)
308  {
309  ASSERTL1(n <= x.num_elements()+x.GetOffset(),"Array out of bounds");
310  ASSERTL1(n <= y.num_elements()+y.GetOffset(),"Array out of bounds");
311 
312  Scatr(n,&x[0],&y[0],&z[0]);
313  }
314 
315 
316  /// \brief Assemble z[y[i]] += x[i]; z should be zero'd first
317  template<class T> void Assmb(int n, const Array<OneD,T> &x, const Array<OneD,int> &y, Array<OneD,T> &z)
318  {
319  ASSERTL1(n <= x.num_elements()+x.GetOffset(),"Array out of bounds");
320  ASSERTL1(n <= y.num_elements()+y.GetOffset(),"Array out of bounds");
321 
322  Assmb(n,&x[0],&y[0],&z[0]);
323  }
324 
325 
326  /************* Reduction routines *****************/
327 
328  /// \brief Subtract return sum(x)
329  template<class T> T Vsum( int n, const Array<OneD, const T> &x, const int incx)
330  {
331  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
332 
333  return Vsum(n,&x[0],incx);
334  }
335 
336 
337  /// \brief Return the index of the maximum element in x
338  template<class T> int Imax( int n, const Array<OneD, const T> &x, const int incx)
339  {
340  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
341 
342  return Imax(n,&x[0],incx);
343  }
344 
345  /// \brief Return the maximum element in x -- called vmax to avoid
346  /// conflict with max
347  template<class T> T Vmax( int n, const Array<OneD, const T> &x, const int incx)
348  {
349  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
350 
351  return Vmax(n,&x[0],incx);
352  }
353 
354  /// \brief Return the index of the maximum absolute element in x
355  template<class T> int Iamax( int n, const Array<OneD, const T> &x, const int incx)
356  {
357  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
358 
359  return Iamax(n,&x[0],incx);
360 
361  }
362 
363  /// \brief Return the maximum absolute element in x
364  /// called vamax to avoid conflict with max
365  template<class T> T Vamax( int n, const Array<OneD, const T> &x, const int incx)
366  {
367  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
368 
369  return Vamax(n,&x[0],incx);
370  }
371 
372 
373  /// \brief Return the index of the minimum element in x
374  template<class T> int Imin( int n, const Array<OneD, const T> &x, const int incx)
375  {
376  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
377 
378  return Imin(n,&x[0],incx);
379  }
380 
381  /// \brief Return the minimum element in x - called vmin to avoid
382  /// conflict with min
383  template<class T> T Vmin( int n, const Array<OneD, const T> &x, const int incx)
384  {
385  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
386 
387  return Vmin(n,&x[0],incx);
388  }
389 
390  /// \brief
391  template<class T> T Dot(int n,
392  const Array<OneD, const T> &w,
393  const Array<OneD, const T> &x)
394  {
395  ASSERTL1(n <= w.num_elements()+w.GetOffset(),"Array out of bounds");
396  ASSERTL1(n <= x.num_elements()+x.GetOffset(),"Array out of bounds");
397 
398  return Dot(n,&w[0],&x[0]);
399  }
400 
401  /// \brief
402  template<class T> T Dot(int n,
403  const Array<OneD, const T> &w, const int incw,
404  const Array<OneD, const T> &x, const int incx)
405  {
406  ASSERTL1(n*incw <= w.num_elements()+w.GetOffset(),"Array out of bounds");
407  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
408 
409  return Dot(n,&w[0],incw,&x[0],incx);
410  }
411 
412  /// \brief
413  template<class T> T Dot2(int n,
414  const Array<OneD, const T> &w,
415  const Array<OneD, const T> &x,
416  const Array<OneD, const int> &y)
417  {
418  ASSERTL1(n <= w.num_elements()+w.GetOffset(),"Array out of bounds");
419  ASSERTL1(n <= x.num_elements()+x.GetOffset(),"Array out of bounds");
420  ASSERTL1(n <= y.num_elements()+y.GetOffset(),"Array out of bounds");
421 
422  return Dot2(n,&w[0],&x[0],&y[0]);
423  }
424 
425  /// \brief
426  template<class T> T Ddot(int n,
427  const Array<OneD, const T> &w, const int incw,
428  const Array<OneD, const T> &x, const int incx,
429  const Array<OneD, const int> &y, const int incy)
430  {
431  ASSERTL1(n*incw <= w.num_elements()+w.GetOffset(),"Array out of bounds");
432  ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Array out of bounds");
433  ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),"Array out of bounds");
434 
435  return Dot2(n,&w[0],incw,&x[0],incx,&y[0],incy);
436  }
437 
438  /********** Memory routines ***********************/
439 
440  template<class T> void Vcopy(int n, const Array<OneD, const T> &x, int incx, Array<OneD,T> &y, int const incy)
441  {
442  ASSERTL1(static_cast<unsigned int>(std::abs(n*incx)) <= x.num_elements()+x.GetOffset(),"Array out of bounds");
443  ASSERTL1(static_cast<unsigned int>(std::abs(n*incy)) <= y.num_elements()+y.GetOffset(),"Array out of bounds");
444 
445  Vcopy(n,&x[0],incx,&y[0],incy);
446  }
447 
448  template<class T> void Reverse(int n, const Array<OneD, const T> &x, int incx, Array<OneD,T> &y, int const incy)
449  {
450  ASSERTL1(static_cast<unsigned int>(std::abs(n*incx)) <= x.num_elements()+x.GetOffset(),"Array out of bounds");
451  ASSERTL1(static_cast<unsigned int>(std::abs(n*incy)) <= y.num_elements()+y.GetOffset(),"Array out of bounds");
452 
453  Reverse(n,&x[0],incx,&y[0],incy);
454  }
455 
456  }
457 #endif //VECTORMATHARRAY_HPP
458 
Generated on Sun Mar 15 2015 00:11:36 for Nektar++ by   doxygen 1.8.1.2