diff options
author | Gael Guennebaud <g.gael@free.fr> | 2010-11-21 10:09:33 +0100 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2010-11-21 10:09:33 +0100 |
commit | 0020ea544a83a8ffb394b8b84bb2a5c865d2dd73 (patch) | |
tree | 8342b505daa006f8fab297ebe220292a28a16d5d /blas | |
parent | 12bfe5e7184a88b5e49c19ad594017119667fb8f (diff) |
implement HEMV level2 blas routine
Diffstat (limited to 'blas')
-rw-r--r-- | blas/level2_impl.h | 37 |
1 files changed, 36 insertions, 1 deletions
diff --git a/blas/level2_impl.h b/blas/level2_impl.h index 92d058a5c..87a4d04df 100644 --- a/blas/level2_impl.h +++ b/blas/level2_impl.h @@ -500,8 +500,43 @@ int EIGEN_BLAS_FUNC(spr2)(char *uplo, int *n, RealScalar *alpha, RealScalar *x, * where alpha and beta are scalars, x and y are n element vectors and * A is an n by n hermitian matrix. */ -int EIGEN_BLAS_FUNC(hemv)(char *uplo, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *x, int *incx, RealScalar *pbeta, RealScalar *y, int *incy) +int EIGEN_BLAS_FUNC(hemv)(char *uplo, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *px, int *incx, RealScalar *pbeta, RealScalar *py, int *incy) { + Scalar* a = reinterpret_cast<Scalar*>(pa); + Scalar* x = reinterpret_cast<Scalar*>(px); + Scalar* y = reinterpret_cast<Scalar*>(py); + Scalar alpha = *reinterpret_cast<Scalar*>(palpha); + Scalar beta = *reinterpret_cast<Scalar*>(pbeta); + + // check arguments + int info = 0; + if(UPLO(*uplo)==INVALID) info = 1; + else if(*n<0) info = 2; + else if(*lda<std::max(1,*n)) info = 5; + else if(*incx==0) info = 7; + else if(*incy==0) info = 10; + if(info) + return xerbla_(SCALAR_SUFFIX_UP"HEMV ",&info,6); + + if(*n==0) + return 1; + + Scalar* actual_x = get_compact_vector(x,*n,*incx); + Scalar* actual_y = get_compact_vector(y,*n,*incy); + + if(beta!=Scalar(1)) + vector(actual_y, *n) *= beta; + + if(alpha!=Scalar(0)) + { + // TODO performs a direct call to the underlying implementation function + if(UPLO(*uplo)==UP) vector(actual_y,*n).noalias() += matrix(a,*n,*n,*lda).selfadjointView<Upper>() * (alpha * vector(actual_x,*n)); + else if(UPLO(*uplo)==LO) vector(actual_y,*n).noalias() += matrix(a,*n,*n,*lda).selfadjointView<Lower>() * (alpha * vector(actual_x,*n)); + } + + if(actual_x!=x) delete[] actual_x; + if(actual_y!=y) delete[] copy_back(actual_y,y,*n,*incy); + return 1; } |