36 #ifndef NEKTAR_LIB_LIBUTILITIES_BASSICUTILS_VECTORMATHARRAY_HPP
37 #define NEKTAR_LIB_LIBUTILITIES_BASSICUTILS_VECTORMATHARRAY_HPP
51 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Out of bounds");
53 Fill(n,alpha,&x[0],incx);
58 const int incx,
int outseed = 9999)
60 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Out of bounds");
68 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
69 ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
70 ASSERTL1(n*incz <= z.num_elements()+z.GetOffset(),
"Array out of bounds");
72 Vmul(n,&x[0],incx,&y[0],incy,&z[0],incz);
77 ASSERTL1(n*incx <= x.num_elements(),
"Array out of bounds");
78 ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
79 ASSERTL1(n*incz <= z.num_elements()+z.GetOffset(),
"Array out of bounds");
81 Vmul(n,x.origin(),incx,&y[0],incy,&z[0],incz);
88 ASSERTL1(static_cast<unsigned int>(n*incx) <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
89 ASSERTL1(static_cast<unsigned int>(n*incy) <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
91 Smul(n,alpha, &x[0],incx,&y[0],incy);
97 ASSERTL1(static_cast<unsigned int>(n*incx) <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
98 ASSERTL1(static_cast<unsigned int>(n*incy) <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
99 ASSERTL1(static_cast<unsigned int>(n*incz) <= z.num_elements()+z.GetOffset(),
"Array out of bounds");
101 Vdiv(n,&x[0],incx,&y[0],incy,&z[0],incz);
108 ASSERTL1(static_cast<unsigned int>(n*incx) <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
109 ASSERTL1(static_cast<unsigned int>(n*incy) <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
111 Sdiv(n,alpha,&x[0],incx,&y[0],incy);
118 ASSERTL1(static_cast<unsigned int>(n*incx) <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
119 ASSERTL1(static_cast<unsigned int>(n*incy) <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
120 ASSERTL1(static_cast<unsigned int>(n*incz) <= z.num_elements()+z.GetOffset(),
"Array out of bounds");
122 Vadd(n,&x[0],incx,&y[0],incy,&z[0],incz);
129 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
130 ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
132 Sadd(n,alpha,&x[0],incx,&y[0],incy);
138 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
139 ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
140 ASSERTL1(n*incz <= z.num_elements()+z.GetOffset(),
"Array out of bounds");
142 Vsub(n,&x[0],incx,&y[0],incy,&z[0],incz);
149 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
158 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
166 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
167 ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
169 Vlog(n, &x[0], incx, &y[0], incy);
175 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
176 ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
178 Vexp(n, &x[0], incx, &y[0], incy);
183 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
184 ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
186 Vpow(n, &x[0], incx, f, &y[0], incy);
192 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
193 ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
195 Vsqrt(n,&x[0],incx,&y[0],incy);
201 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
202 ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
204 Vabs(n,&x[0],incx,&y[0],incy);
210 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)
212 ASSERTL1(n*incw <= w.num_elements()+w.GetOffset(),
"Array out of bounds");
213 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
214 ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
215 ASSERTL1(n*incz <= z.num_elements()+z.GetOffset(),
"Array out of bounds");
217 Vvtvp(n,&w[0],incw,&x[0],incx,&y[0],incy,&z[0],incz);
220 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)
222 ASSERTL1(n*incw <= w.num_elements(),
"Array out of bounds");
223 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
224 ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
225 ASSERTL1(n*incz <= z.num_elements()+z.GetOffset(),
"Array out of bounds");
227 Vvtvp(n,w.origin(),incw,&x[0],incx,&y[0],incy,&z[0],incz);
231 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)
233 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
234 ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
235 ASSERTL1(n*incz <= z.num_elements()+z.GetOffset(),
"Array out of bounds");
237 Svtvp(n,alpha,&x[0],incx,&y[0],incy,&z[0],incz);
242 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)
244 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
245 ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
246 ASSERTL1(n*incz <= z.num_elements()+z.GetOffset(),
"Array out of bounds");
248 Svtvm(n,alpha,&x[0],incx,&y[0],incy,&z[0],incz);
253 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)
255 ASSERTL1(n*incw <= w.num_elements()+w.GetOffset(),
"Array out of bounds");
256 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
257 ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
258 ASSERTL1(n*incz <= z.num_elements()+z.GetOffset(),
"Array out of bounds");
260 Vvtvm(n,&w[0],incw,&x[0],incx,&y[0],incy,&z[0],incz);
273 ASSERTL1(n*incv <= v.num_elements()+v.GetOffset(),
"Array out of bounds");
274 ASSERTL1(n*incw <= w.num_elements()+w.GetOffset(),
"Array out of bounds");
275 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
276 ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
277 ASSERTL1(n*incz <= z.num_elements()+z.GetOffset(),
"Array out of bounds");
279 Vvtvvtp(n,&v[0],incv,&w[0],incw,&x[0],incx,&y[0],incy,&z[0],incz);
283 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)
285 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
286 ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
287 ASSERTL1(n*incz <= z.num_elements()+z.GetOffset(),
"Array out of bounds");
289 Svtsvtp(n,alpha,&x[0],incx,beta,&y[0],incy,&z[0],incz);
301 ASSERTL1(n <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
302 ASSERTL1(n <= z.num_elements()+z.GetOffset(),
"Array out of bounds");
304 Gathr(n,&x[0],&y[0],&z[0]);
311 ASSERTL1(n <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
312 ASSERTL1(n <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
314 Scatr(n,&x[0],&y[0],&z[0]);
321 ASSERTL1(n <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
322 ASSERTL1(n <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
324 Assmb(n,&x[0],&y[0],&z[0]);
333 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
335 return Vsum(n,&x[0],incx);
342 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
344 return Imax(n,&x[0],incx);
351 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
353 return Vmax(n,&x[0],incx);
359 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
361 return Iamax(n,&x[0],incx);
369 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
371 return Vamax(n,&x[0],incx);
378 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
380 return Imin(n,&x[0],incx);
387 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
389 return Vmin(n,&x[0],incx);
395 ASSERTL1(n * incx <= x.num_elements() + x.GetOffset(),
"Array out of bounds");
397 return Nnan(n, &x[0], incx);
401 template<
class T> T
Dot(
int n,
405 ASSERTL1(n <= w.num_elements()+w.GetOffset(),
"Array out of bounds");
406 ASSERTL1(n <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
408 return Dot(n,&w[0],&x[0]);
412 template<
class T> T
Dot(
int n,
416 ASSERTL1(n*incw <= w.num_elements()+w.GetOffset(),
"Array out of bounds");
417 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
419 return Dot(n,&w[0],incw,&x[0],incx);
423 template<
class T> T
Dot2(
int n,
428 ASSERTL1(n <= w.num_elements()+w.GetOffset(),
"Array out of bounds");
429 ASSERTL1(n <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
430 ASSERTL1(n <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
432 return Dot2(n,&w[0],&x[0],&y[0]);
436 template<
class T> T
Ddot(
int n,
441 ASSERTL1(n*incw <= w.num_elements()+w.GetOffset(),
"Array out of bounds");
442 ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
443 ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
445 return Dot2(n,&w[0],incw,&x[0],incx,&y[0],incy);
452 ASSERTL1(static_cast<unsigned int>(std::abs(n*incx)) <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
453 ASSERTL1(static_cast<unsigned int>(std::abs(n*incy)) <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
455 Vcopy(n,&x[0],incx,&y[0],incy);
460 ASSERTL1(static_cast<unsigned int>(std::abs(n*incx)) <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
461 ASSERTL1(static_cast<unsigned int>(std::abs(n*incy)) <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
463 Reverse(n,&x[0],incx,&y[0],incy);
467 #endif //VECTORMATHARRAY_HPP
void Vpow(int n, const T *x, const int incx, const T f, T *y, const int incy)
void Gathr(int n, const T *x, const int *y, T *z)
Gather vector z[i] = x[y[i]].
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
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
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
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/y.
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.
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
int Imax(int n, const T *x, const int incx)
Return the index of the maximum element in x.
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
int Nnan(int n, const T *x, const int incx)
Return number of NaN elements of x.
void Vexp(int n, const T *x, const int incx, T *y, const int incy)
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
void Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
void Neg(int n, T *x, const int incx)
Negate x = -x.
void Assmb(int n, const T *x, const int *y, T *z)
Assemble z[y[i]] += x[i]; z should be zero'd first.
T Vamax(int n, const T *x, const int incx)
Return the maximum absolute element in x called vamax to avoid conflict with max. ...
int Iamax(int n, const T *x, const int incx)
Return the index of the maximum absolute element in x.
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
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)
T Dot(int n, const T *w, const T *x)
vvtvp (vector times vector times vector): z = w*x*y
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.
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):
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
T Dot2(int n, const T *w, const T *x, const int *y)
vvtvp (vector times vector times vector): z = w*x*y
void FillWhiteNoise(int n, const T eps, T *x, const int incx, int outseed)
Fills a vector with white noise.
void Vlog(int n, const T *x, const int incx, T *y, const int incy)
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):
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
void Zero(int n, T *x, const int incx)
Zero vector.
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
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.
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.