From 68f8519327edc27003e7ca69616759517cdbbbec Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Fri, 19 Nov 2010 16:51:52 +0100 Subject: implement HER and HER2 blas routines --- blas/common.h | 12 +++++++ blas/level2_impl.h | 91 +++++++++++++++++++++++++++++++++++++----------------- 2 files changed, 74 insertions(+), 29 deletions(-) (limited to 'blas') diff --git a/blas/common.h b/blas/common.h index d56815ce3..9eaca2c2f 100644 --- a/blas/common.h +++ b/blas/common.h @@ -128,6 +128,18 @@ Map > vector(T* data, int size) return Map >(data, size); } +template +T* get_compact_vector(T* x, int n, int incx) +{ + if(incx==1) + return x; + + T* ret = new Scalar[n]; + if(incx<0) vector(ret,n) = vector(x,n,-incx).reverse(); + else vector(ret,n) = vector(x,n, incx); + return ret; +} + #define EIGEN_BLAS_FUNC(X) EIGEN_CAT(SCALAR_SUFFIX,X##_) #endif // EIGEN_BLAS_COMMON_H diff --git a/blas/level2_impl.h b/blas/level2_impl.h index 2dc059c14..7eaf83835 100644 --- a/blas/level2_impl.h +++ b/blas/level2_impl.h @@ -239,20 +239,13 @@ int EIGEN_BLAS_FUNC(syr)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, return 1; // if the increment is not 1, let's copy it to a temporary vector to enable vectorization - Scalar* x_cpy = x; - if(*incx!=1) - { - x_cpy = new Scalar[*n]; - if(*incx<0) vector(x_cpy,*n) = vector(x,*n,-*incx).reverse(); - else vector(x_cpy,*n) = vector(x,*n,*incx); - } + Scalar* x_cpy = get_compact_vector(x,*n,*incx); // TODO perform direct calls to underlying implementation if(UPLO(*uplo)==LO) matrix(c,*n,*n,*ldc).selfadjointView().rankUpdate(vector(x_cpy,*n), alpha); else if(UPLO(*uplo)==UP) matrix(c,*n,*n,*ldc).selfadjointView().rankUpdate(vector(x_cpy,*n), alpha); - if(*incx!=1) - delete[] x_cpy; + if(x_cpy!=x) delete[] x_cpy; // func[code](*n, a, *inca, c, *ldc, alpha); return 1; @@ -294,29 +287,15 @@ int EIGEN_BLAS_FUNC(syr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px if(alpha==Scalar(0)) return 1; - // if the increment is not 1, let's copy it to a temporary vector to enable vectorization - Scalar* x_cpy = x; - if(*incx!=1) - { - x_cpy = new Scalar[*n]; - if(*incx<0) vector(x_cpy,*n) = vector(x,*n,-*incx).reverse(); - else vector(x_cpy,*n) = vector(x,*n, *incx); - } - - Scalar* y_cpy = y; - if(*incy!=1) - { - y_cpy = new Scalar[*n]; - if(*incy<0) vector(y_cpy,*n) = vector(y,*n,-*incy).reverse(); - else vector(y_cpy,*n) = vector(y,*n, *incy); - } + Scalar* x_cpy = get_compact_vector(x,*n,*incx); + Scalar* y_cpy = get_compact_vector(y,*n,*incy); // TODO perform direct calls to underlying implementation if(UPLO(*uplo)==LO) matrix(c,*n,*n,*ldc).selfadjointView().rankUpdate(vector(x_cpy,*n), vector(y_cpy,*n), alpha); else if(UPLO(*uplo)==UP) matrix(c,*n,*n,*ldc).selfadjointView().rankUpdate(vector(x_cpy,*n), vector(y_cpy,*n), alpha); - if(*incx!=1) delete[] x_cpy; - if(*incy!=1) delete[] y_cpy; + if(x_cpy!=x) delete[] x_cpy; + if(y_cpy!=y) delete[] y_cpy; // int code = UPLO(*uplo); // if(code>=2 || func[code]==0) @@ -523,8 +502,33 @@ int EIGEN_BLAS_FUNC(hpr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *x, * where alpha is a real scalar, x is an n element vector and A is an * n by n hermitian matrix. */ -int EIGEN_BLAS_FUNC(her)(char *uplo, int *n, RealScalar *alpha, RealScalar *x, int *incx, RealScalar *a, int *lda) +int EIGEN_BLAS_FUNC(her)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *pa, int *lda) { + Scalar* x = reinterpret_cast(px); + Scalar* a = reinterpret_cast(pa); + RealScalar alpha = *reinterpret_cast(palpha); + + int info = 0; + if(UPLO(*uplo)==INVALID) info = 1; + else if(*n<0) info = 2; + else if(*incx==0) info = 5; + else if(*lda().rankUpdate(vector(x_cpy,*n), alpha); + else if(UPLO(*uplo)==UP) matrix(a,*n,*n,*lda).selfadjointView().rankUpdate(vector(x_cpy,*n), alpha); + + matrix(a,*n,*n,*lda).diagonal().imag().setZero(); + + if(x_cpy!=x) delete[] x_cpy; + return 1; } @@ -535,8 +539,37 @@ int EIGEN_BLAS_FUNC(her)(char *uplo, int *n, RealScalar *alpha, RealScalar *x, i * where alpha is a scalar, x and y are n element vectors and A is an n * by n hermitian matrix. */ -int EIGEN_BLAS_FUNC(her2)(char *uplo, int *n, RealScalar *alpha, RealScalar *x, int *incx, RealScalar *y, int *incy, RealScalar *a, int *lda) +int EIGEN_BLAS_FUNC(her2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pa, int *lda) { + Scalar* x = reinterpret_cast(px); + Scalar* y = reinterpret_cast(py); + Scalar* a = reinterpret_cast(pa); + Scalar alpha = *reinterpret_cast(palpha); + + int info = 0; + if(UPLO(*uplo)==INVALID) info = 1; + else if(*n<0) info = 2; + else if(*incx==0) info = 5; + else if(*incy==0) info = 7; + else if(*lda().rankUpdate(vector(x_cpy,*n),vector(y_cpy,*n),alpha); + else if(UPLO(*uplo)==UP) matrix(a,*n,*n,*lda).selfadjointView().rankUpdate(vector(x_cpy,*n),vector(y_cpy,*n),alpha); + + matrix(a,*n,*n,*lda).diagonal().imag().setZero(); + + if(x_cpy!=x) delete[] x_cpy; + if(y_cpy!=y) delete[] y_cpy; + return 1; } -- cgit v1.2.3