45 template<
class T>
void Fill(
int n,
const T alpha, T *x,
const int incx )
58 #define IM1 2147483563
59 #define IM2 2147483399
69 #define NDIV (1+IMM1/NTAB)
71 #define RNMX (1.0-EPS)
74 template<
class T> T
ran2 (
long* idum)
85 static long idum2=123456789;
91 if (-(*idum) < 1) *idum = 1;
92 else *idum = -(*idum);
94 for (j=
NTAB+7; j>=0; j--) {
97 if (*idum < 0) *idum +=
IM1;
98 if (j <
NTAB) iv[j] = *idum;
105 if (*idum < 0) *idum +=
IM1;
109 if (idum2 < 0) idum2 +=
IM2;
114 if (iy < 1) iy +=
IMM1;
135 #ifdef NEKTAR_USE_THREAD_SAFETY
136 static boost::mutex mutex;
143 const int incx,
int outseed)
145 #ifdef NEKTAR_USE_THREAD_SAFETY
147 boost::mutex::scoped_lock l(mutex);
153 static long seed = 0;
158 seed = long(outseed);
169 v1 = 2.0 * ran2<T> (&seed) - 1.0;
170 v2 = 2.0 * ran2<T> (&seed) - 1.0;
172 }
while (rsq >= 1.0 || rsq == 0.0);
173 fac =
sqrt(-2.0 *
log (rsq) / rsq);
189 const int incx,
int outseed);
192 template<
class T>
void Vmul(
int n,
const T *x,
const int incx,
const T *y,
193 const int incy, T *z,
const int incz)
196 if (incx == 1 && incy == 1 && incz == 1)
225 template<
class T>
void Smul(
int n,
const T alpha,
const T *x,
const int incx,
226 T *y,
const int incy)
229 if (incx == 1 && incy == 1)
257 template<
class T>
void Vdiv(
int n,
const T *x,
const int incx,
const T *y,
258 const int incy, T*z,
const int incz)
261 if (incx == 1 && incy == 1)
291 template<
class T>
void Sdiv(
int n,
const T alpha,
const T *x,
292 const int incx, T *y,
const int incy)
295 if (incx == 1 && incy == 1)
322 template<
class T>
void Vadd(
int n,
const T *x,
const int incx,
const T *y,
323 const int incy, T *z,
const int incz)
341 template<
class T>
void Sadd(
int n,
const T alpha,
const T *x,
342 const int incx, T *y,
const int incy)
345 if (incx == 1 && incy == 1)
372 template<
class T>
void Vsub(
int n,
const T *x,
const int incx,
const T *y,
373 const int incy, T *z,
const int incz)
376 if (incx == 1 && incy == 1 && incz == 1)
405 template<
class T>
void Ssub(
int n,
const T alpha,
const T *x,
406 const int incx, T *y,
const int incy)
409 if (incx == 1 && incy == 1)
436 template<
class T>
void Zero(
int n, T *x,
const int incx)
440 std::memset(x,
'\0', n*
sizeof(T));
461 template<
class T>
void Neg(
int n, T *x,
const int incx)
475 template<
class T>
void Vsqrt(
int n,
const T *x,
const int incx,
476 T *y,
const int incy)
493 template<
class T>
void Vabs(
int n,
const T *x,
const int incx,
494 T *y,
const int incy)
498 *y = ( *x >0)? *x:-(*x);
514 const T *w,
const int incw,
515 const T *x,
const int incx,
516 const T *y,
const int incy,
517 T *z,
const int incz)
521 *z = (*w) * (*x) + (*y);
541 template<
class T>
void Vvtvm(
int n,
const T *w,
const int incw,
const T *x,
542 const int incx,
const T *y,
const int incy,
543 T *z,
const int incz)
547 *z = (*w) * (*x) - (*y);
566 const int incx,
const T *y,
const int incy,
567 T *z,
const int incz)
570 if (incx == 1 && incy == 1 && incz == 1)
574 *z = alpha * (*x) + (*y);
584 *z = alpha * (*x) + (*y);
602 template<
class T>
void Svtvm(
int n,
const T alpha,
const T *x,
603 const int incx,
const T *y,
const int incy,
604 T *z,
const int incz)
608 *z = alpha * (*x) - (*y);
626 const T* v,
int incv,
627 const T* w,
int incw,
628 const T* x,
int incx,
629 const T* y,
int incy,
634 *z = (*v) * (*w) + (*x) * (*y);
659 const T* v,
int incv,
660 const T* w,
int incw,
661 const T* x,
int incx,
662 const T* y,
int incy,
667 *z = (*v) * (*w) - (*x) * (*y);
693 const T* x,
int incx,
695 const T* y,
int incy,
700 *z = alpha * (*x) + beta * (*y);
725 const T* v,
int incv,
726 const T* w,
int incw,
727 const T* x,
int incx,
732 *z = alpha * (*v) + (*w) + (*x);
756 template<
class T>
void Gathr(
int n,
const T *
sign,
const T *x,
const int *y,
761 *z++ = *(
sign++) * (*(x + *y++));
772 template<
class T>
void Scatr(
int n,
const T *x,
const int *y,
777 *(z + *(y++)) = *(x++);
787 template<
class T>
void Scatr(
int n,
const T *
sign,
const T *x,
const int *y,
794 *(z + *(y++)) = *(
sign++) * (*(x++));
813 template<
class T>
void Assmb(
int n,
const T *x,
const int *y,
818 *(z + *(y++)) += *(x++);
828 template<
class T>
void Assmb(
int n,
const T *
sign,
const T *x,
const int *y,
833 *(z + *(y++)) += *(
sign++) * (*(x++));
846 template<
class T> T
Vsum(
int n,
const T *x,
const int incx)
866 template<
class T>
int Imax(
int n,
const T *x,
const int incx)
869 int i, indx = ( n > 0 ) ? 0 : -1;
872 for (i = 0; i < n; i++)
892 template<
class T> T
Vmax(
int n,
const T *x,
const int incx)
915 template<
class T>
int Iamax(
int n,
const T *x,
const int incx)
918 int i, indx = ( n > 0 ) ? 0 : -1;
922 for (i = 0; i < n; i++)
924 xm = (*x > 0)? *x: -*x;
942 template<
class T> T
Vamax(
int n,
const T *x,
const int incx)
950 xm = (*x > 0)? *x: -*x;
966 template<
class T>
int Imin(
int n,
const T *x,
const int incx)
969 int i, indx = ( n > 0 ) ? 0 : -1;
992 template<
class T> T
Vmin(
int n,
const T *x,
const int incx)
1015 template<
class T>
int Nnan(
int n,
const T *x,
const int incx)
1038 template<
class T> T
Dot(
int n,
1061 template<
class T> T
Dot(
int n,
1062 const T *w,
const int incw,
1063 const T *x,
const int incx)
1093 sum += (*y == 1 ? (*w) * (*x) : 0 );
1113 const T *w,
const int incw,
1114 const T *x,
const int incx,
1115 const int *y,
const int incy)
1121 sum += (*y == 1 ? (*w) * (*x) : 0.0 );
1132 const int *y,
const int incy);
1136 const int *y,
const int incy);
1198 template<
typename T>
1199 void Vcopy(
int n,
const T *x,
const int incx,
1200 T *y,
const int incy)
1202 if( incx ==1 && incy == 1)
1204 memcpy(y,x,n*
sizeof(T));
1226 template<
class T>
void Reverse(
int n,
const T *x,
const int incx, T *y,
const int incy)
1236 y[nloop] = x[nloop];
1238 const T* x_end = x + (n-1)*incx;
1239 T* y_end = y + (n-1)*incy;
1240 for (i = 0; i < nloop; ++i) {
#define LIB_UTILITIES_EXPORT
#define sign(a, b)
return the sign(b)*a
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
void Ssub(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha - x.
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):
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.
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 Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
T Dot2(int n, const T *w, const T *x, const int *y)
vvtvp (vector times vector times vector): z = w*x*y
void Neg(int n, T *x, const int incx)
Negate x = -x.
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
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
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
void Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
T Dot(int n, const T *w, const T *x)
vvtvp (vector times vector times vector): z = w*x*y
void Gathr(int n, const T *sign, const T *x, const int *y, T *z)
Gather vector z[i] = sign[i]*x[y[i]].
void Assmb(int n, const T *x, const int *y, T *z)
Assemble z[y[i]] += x[i]; z should be zero'd first.
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 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 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
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):
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/y.
int Imax(int n, const T *x, const int incx)
Return the index of the maximum element in x.
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.
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
void Zero(int n, T *x, const int incx)
Zero vector.
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
T ran2(long *idum)
Generates a number from ~Normal(0,1)
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):
void FillWhiteNoise(int n, const T eps, T *x, const int incx, int outseed)
Fills a vector with white noise.
int Nnan(int n, const T *x, const int incx)
Return number of NaN elements of x.
T Vamax(int n, const T *x, const int incx)
Return the maximum absolute element in x called vamax to avoid conflict with max.
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha - x.
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
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 Vcopy(int n, const T *x, const int incx, T *y, const int incy)
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.
int Iamax(int n, const T *x, const int incx)
Return the index of the maximum absolute element in x.
scalarT< T > log(scalarT< T > in)
scalarT< T > sqrt(scalarT< T > in)