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             ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Out of bounds");
 
   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");
 
   70             Vmul(n,&x[0],incx,&y[0],incy,&z[0],incz);
 
   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");
 
   79             Vmul(n,x.origin(),incx,&y[0],incy,&z[0],incz);
 
   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");
 
   89             Smul(n,alpha, &x[0],incx,&y[0],incy);
 
   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");
 
   99             Vdiv(n,&x[0],incx,&y[0],incy,&z[0],incz);
 
  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");
 
  109             Sdiv(n,alpha,&x[0],incx,&y[0],incy);
 
  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");
 
  120             Vadd(n,&x[0],incx,&y[0],incy,&z[0],incz);
 
  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");
 
  130             Sadd(n,alpha,&x[0],incx,&y[0],incy);
 
  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");
 
  140             Vsub(n,&x[0],incx,&y[0],incy,&z[0],incz);
 
  147             ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
 
  156             ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
 
  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");
 
  167             Vlog(n, &x[0], incx, &y[0], incy);
 
  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");
 
  176             Vexp(n, &x[0], incx, &y[0], incy);
 
  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");
 
  184             Vpow(n, &x[0], incx, f, &y[0], incy);
 
  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");
 
  193             Vsqrt(n,&x[0],incx,&y[0],incy);
 
  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");
 
  202             Vabs(n,&x[0],incx,&y[0],incy);
 
  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)
 
  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");
 
  215             Vvtvp(n,&w[0],incw,&x[0],incx,&y[0],incy,&z[0],incz);
 
  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)
 
  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");
 
  225             Vvtvp(n,w.origin(),incw,&x[0],incx,&y[0],incy,&z[0],incz);
 
  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)
 
  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");
 
  235             Svtvp(n,alpha,&x[0],incx,&y[0],incy,&z[0],incz);
 
  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)
 
  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");
 
  246             Svtvm(n,alpha,&x[0],incx,&y[0],incy,&z[0],incz);
 
  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)
 
  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");
 
  258             Vvtvm(n,&w[0],incw,&x[0],incx,&y[0],incy,&z[0],incz);
 
  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");
 
  277             Vvtvvtp(n,&v[0],incv,&w[0],incw,&x[0],incx,&y[0],incy,&z[0],incz);
 
  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)
 
  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");
 
  287             Svtsvtp(n,alpha,&x[0],incx,beta,&y[0],incy,&z[0],incz);
 
  299             ASSERTL1(n <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
 
  300             ASSERTL1(n <= z.num_elements()+z.GetOffset(),
"Array out of bounds");
 
  302             Gathr(n,&x[0],&y[0],&z[0]);
 
  309             ASSERTL1(n <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
 
  310             ASSERTL1(n <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
 
  312             Scatr(n,&x[0],&y[0],&z[0]);
 
  319             ASSERTL1(n <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
 
  320             ASSERTL1(n <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
 
  322             Assmb(n,&x[0],&y[0],&z[0]);
 
  331             ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
 
  333             return Vsum(n,&x[0],incx);
 
  340             ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
 
  342             return Imax(n,&x[0],incx);
 
  349             ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
 
  351             return Vmax(n,&x[0],incx);
 
  357             ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
 
  359             return Iamax(n,&x[0],incx);
 
  367             ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
 
  369             return Vamax(n,&x[0],incx);            
 
  376             ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
 
  378             return Imin(n,&x[0],incx);
 
  385             ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
 
  387             return Vmin(n,&x[0],incx);
 
  393             ASSERTL1(n * incx <= x.num_elements() + x.GetOffset(), 
"Array out of bounds");
 
  395             return Nnan(n, &x[0], incx);
 
  399         template<
class T> T 
Dot(
int n,
 
  403             ASSERTL1(n <= w.num_elements()+w.GetOffset(),
"Array out of bounds");
 
  404             ASSERTL1(n <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
 
  406             return Dot(n,&w[0],&x[0]);
 
  410         template<
class T> T 
Dot(
int n,
 
  414             ASSERTL1(n*incw <= w.num_elements()+w.GetOffset(),
"Array out of bounds");
 
  415             ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
 
  417             return Dot(n,&w[0],incw,&x[0],incx);
 
  421         template<
class T> T 
Dot2(
int n,
 
  426             ASSERTL1(n <= w.num_elements()+w.GetOffset(),
"Array out of bounds");
 
  427             ASSERTL1(n <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
 
  428             ASSERTL1(n <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
 
  430             return Dot2(n,&w[0],&x[0],&y[0]);
 
  434         template<
class T> T 
Ddot(
int n,
 
  439             ASSERTL1(n*incw <= w.num_elements()+w.GetOffset(),
"Array out of bounds");
 
  440             ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
 
  441             ASSERTL1(n*incy <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
 
  443             return Dot2(n,&w[0],incw,&x[0],incx,&y[0],incy);
 
  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");
 
  453             Vcopy(n,&x[0],incx,&y[0],incy);
 
  458             ASSERTL1(static_cast<unsigned int>(std::abs(n*incx)) <= x.num_elements()+x.GetOffset(),
"Array out of bounds");
 
  459             ASSERTL1(static_cast<unsigned int>(std::abs(n*incy)) <= y.num_elements()+y.GetOffset(),
"Array out of bounds");
 
  461             Reverse(n,&x[0],incx,&y[0],incy);
 
  465 #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.