From f5f288b741b173a271b9c939ac5231639135dd93 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 22 Nov 2010 18:49:12 +0100 Subject: split level 1 and 2 implementation files into smaller ones and fix a couple of numerical and tricky issues discovered by the lapack test suite --- blas/level2_impl.h | 505 ++++------------------------------------------------- 1 file changed, 36 insertions(+), 469 deletions(-) (limited to 'blas/level2_impl.h') diff --git a/blas/level2_impl.h b/blas/level2_impl.h index 87a4d04df..0e4649f57 100644 --- a/blas/level2_impl.h +++ b/blas/level2_impl.h @@ -41,7 +41,7 @@ int EIGEN_BLAS_FUNC(gemv)(char *opa, int *m, int *n, RealScalar *palpha, RealSca init = true; } - + Scalar* a = reinterpret_cast(pa); Scalar* b = reinterpret_cast(pb); Scalar* c = reinterpret_cast(pc); @@ -59,7 +59,7 @@ int EIGEN_BLAS_FUNC(gemv)(char *opa, int *m, int *n, RealScalar *palpha, RealSca if(info) return xerbla_(SCALAR_SUFFIX_UP"GEMV ",&info,6); - if(*m==0 || *n==0) + if(*m==0 || *n==0 || (alpha==Scalar(0) && beta==Scalar(1))) return 0; int actual_m = *m; @@ -69,10 +69,13 @@ int EIGEN_BLAS_FUNC(gemv)(char *opa, int *m, int *n, RealScalar *palpha, RealSca Scalar* actual_b = get_compact_vector(b,actual_n,*incb); Scalar* actual_c = get_compact_vector(c,actual_m,*incc); - + if(beta!=Scalar(1)) - vector(actual_c, actual_m) *= beta; - + { + if(beta==Scalar(0)) vector(actual_c, actual_m).setZero(); + else vector(actual_c, actual_m) *= beta; + } + int code = OP(*opa); func[code](actual_m, actual_n, a, *lda, actual_b, 1, actual_c, 1, alpha); @@ -131,7 +134,7 @@ int EIGEN_BLAS_FUNC(trsv)(char *uplo, char *opa, char *diag, int *n, RealScalar func[code](*n, a, *lda, actual_b); if(actual_b!=b) delete[] copy_back(actual_b,b,*n,*incb); - + return 0; } @@ -195,147 +198,8 @@ int EIGEN_BLAS_FUNC(trmv)(char *uplo, char *opa, char *diag, int *n, RealScalar copy_back(res.data(),b,*n,*incb); if(actual_b!=b) delete[] actual_b; - - return 0; -} - -// y = alpha*A*x + beta*y -int EIGEN_BLAS_FUNC(symv) (char *uplo, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *px, int *incx, RealScalar *pbeta, RealScalar *py, int *incy) -{ - Scalar* a = reinterpret_cast(pa); - Scalar* x = reinterpret_cast(px); - Scalar* y = reinterpret_cast(py); - Scalar alpha = *reinterpret_cast(palpha); - Scalar beta = *reinterpret_cast(pbeta); - - // check arguments - int info = 0; - if(UPLO(*uplo)==INVALID) info = 1; - else if(*n<0) info = 2; - else if(*lda() * (alpha * vector(actual_x,*n)); - else if(UPLO(*uplo)==LO) vector(actual_y,*n).noalias() += matrix(a,*n,*n,*lda).selfadjointView() * (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; -} - -// C := alpha*x*x' + C -int EIGEN_BLAS_FUNC(syr)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *pc, int *ldc) -{ - -// typedef void (*functype)(int, const Scalar *, int, Scalar *, int, Scalar); -// static functype func[2]; - -// static bool init = false; -// if(!init) -// { -// for(int k=0; k<2; ++k) -// func[k] = 0; -// -// func[UP] = (internal::selfadjoint_product::run); -// func[LO] = (internal::selfadjoint_product::run); - -// init = true; -// } - - Scalar* x = reinterpret_cast(px); - Scalar* c = reinterpret_cast(pc); - 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(*ldc().rankUpdate(vector(x_cpy,*n), alpha); - else if(UPLO(*uplo)==UP) matrix(c,*n,*n,*ldc).selfadjointView().rankUpdate(vector(x_cpy,*n), alpha); - - if(x_cpy!=x) delete[] x_cpy; - -// func[code](*n, a, *inca, c, *ldc, alpha); - return 1; -} - - -// C := alpha*x*y' + alpha*y*x' + C -int EIGEN_BLAS_FUNC(syr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pc, int *ldc) -{ -// typedef void (*functype)(int, const Scalar *, int, const Scalar *, int, Scalar *, int, Scalar); -// static functype func[2]; -// -// static bool init = false; -// if(!init) -// { -// for(int k=0; k<2; ++k) -// func[k] = 0; -// -// func[UP] = (internal::selfadjoint_product::run); -// func[LO] = (internal::selfadjoint_product::run); -// -// init = true; -// } - - Scalar* x = reinterpret_cast(px); - Scalar* y = reinterpret_cast(py); - Scalar* c = reinterpret_cast(pc); - 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(*ldc().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(x_cpy!=x) delete[] x_cpy; - if(y_cpy!=y) delete[] y_cpy; - -// int code = UPLO(*uplo); -// if(code>=2 || func[code]==0) -// return 0; -// func[code](*n, a, *inca, b, *incb, c, *ldc, alpha); - return 1; + return 0; } /** DGBMV performs one of the matrix-vector operations @@ -345,23 +209,12 @@ int EIGEN_BLAS_FUNC(syr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px * where alpha and beta are scalars, x and y are vectors and A is an * m by n band matrix, with kl sub-diagonals and ku super-diagonals. */ -int EIGEN_BLAS_FUNC(gbmv)(char *trans, int *m, int *n, int *kl, int *ku, RealScalar *alpha, RealScalar *a, int *lda, - RealScalar *x, int *incx, RealScalar *beta, RealScalar *y, int *incy) -{ - return 1; -} -/** DSBMV performs the matrix-vector operation - * - * y := alpha*A*x + beta*y, - * - * where alpha and beta are scalars, x and y are n element vectors and - * A is an n by n symmetric band matrix, with k super-diagonals. - */ -int EIGEN_BLAS_FUNC(sbmv)( char *uplo, int *n, int *k, RealScalar *alpha, RealScalar *a, int *lda, - RealScalar *x, int *incx, RealScalar *beta, RealScalar *y, int *incy) -{ - return 1; -} +// int EIGEN_BLAS_FUNC(gbmv)(char *trans, int *m, int *n, int *kl, int *ku, RealScalar *alpha, RealScalar *a, int *lda, +// RealScalar *x, int *incx, RealScalar *beta, RealScalar *y, int *incy) +// { +// return 1; +// } + /** DTBMV performs one of the matrix-vector operations * @@ -370,10 +223,10 @@ int EIGEN_BLAS_FUNC(sbmv)( char *uplo, int *n, int *k, RealScalar *alpha, RealSc * where x is an n element vector and A is an n by n unit, or non-unit, * upper or lower triangular band matrix, with ( k + 1 ) diagonals. */ -int EIGEN_BLAS_FUNC(tbmv)(char *uplo, char *trans, char *diag, int *n, int *k, RealScalar *a, int *lda, RealScalar *x, int *incx) -{ - return 1; -} +// int EIGEN_BLAS_FUNC(tbmv)(char *uplo, char *trans, char *diag, int *n, int *k, RealScalar *a, int *lda, RealScalar *x, int *incx) +// { +// return 1; +// } /** DTBSV solves one of the systems of equations * @@ -386,23 +239,10 @@ int EIGEN_BLAS_FUNC(tbmv)(char *uplo, char *trans, char *diag, int *n, int *k, R * No test for singularity or near-singularity is included in this * routine. Such tests must be performed before calling this routine. */ -int EIGEN_BLAS_FUNC(tbsv)(char *uplo, char *trans, char *diag, int *n, int *k, RealScalar *a, int *lda, RealScalar *x, int *incx) -{ - return 1; -} - -/** DSPMV performs the matrix-vector operation - * - * y := alpha*A*x + beta*y, - * - * where alpha and beta are scalars, x and y are n element vectors and - * A is an n by n symmetric matrix, supplied in packed form. - * - */ -int EIGEN_BLAS_FUNC(spmv)(char *uplo, int *n, RealScalar *alpha, RealScalar *ap, RealScalar *x, int *incx, RealScalar *beta, RealScalar *y, int *incy) -{ - return 1; -} +// int EIGEN_BLAS_FUNC(tbsv)(char *uplo, char *trans, char *diag, int *n, int *k, RealScalar *a, int *lda, RealScalar *x, int *incx) +// { +// return 1; +// } /** DTPMV performs one of the matrix-vector operations * @@ -411,10 +251,10 @@ int EIGEN_BLAS_FUNC(spmv)(char *uplo, int *n, RealScalar *alpha, RealScalar *ap, * where x is an n element vector and A is an n by n unit, or non-unit, * upper or lower triangular matrix, supplied in packed form. */ -int EIGEN_BLAS_FUNC(tpmv)(char *uplo, char *trans, char *diag, int *n, RealScalar *ap, RealScalar *x, int *incx) -{ - return 1; -} +// int EIGEN_BLAS_FUNC(tpmv)(char *uplo, char *trans, char *diag, int *n, RealScalar *ap, RealScalar *x, int *incx) +// { +// return 1; +// } /** DTPSV solves one of the systems of equations * @@ -426,10 +266,10 @@ int EIGEN_BLAS_FUNC(tpmv)(char *uplo, char *trans, char *diag, int *n, RealScala * No test for singularity or near-singularity is included in this * routine. Such tests must be performed before calling this routine. */ -int EIGEN_BLAS_FUNC(tpsv)(char *uplo, char *trans, char *diag, int *n, RealScalar *ap, RealScalar *x, int *incx) -{ - return 1; -} +// int EIGEN_BLAS_FUNC(tpsv)(char *uplo, char *trans, char *diag, int *n, RealScalar *ap, RealScalar *x, int *incx) +// { +// return 1; +// } /** DGER performs the rank 1 operation * @@ -444,7 +284,7 @@ int EIGEN_BLAS_FUNC(ger)(int *m, int *n, Scalar *palpha, Scalar *px, int *incx, Scalar* y = reinterpret_cast(py); Scalar* a = reinterpret_cast(pa); Scalar alpha = *reinterpret_cast(palpha); - + int info = 0; if(*m<0) info = 1; else if(*n<0) info = 2; @@ -453,293 +293,20 @@ int EIGEN_BLAS_FUNC(ger)(int *m, int *n, Scalar *palpha, Scalar *px, int *incx, else if(*lda(pa); - Scalar* x = reinterpret_cast(px); - Scalar* y = reinterpret_cast(py); - Scalar alpha = *reinterpret_cast(palpha); - Scalar beta = *reinterpret_cast(pbeta); - // check arguments - int info = 0; - if(UPLO(*uplo)==INVALID) info = 1; - else if(*n<0) info = 2; - else if(*lda() * (alpha * vector(actual_x,*n)); - else if(UPLO(*uplo)==LO) vector(actual_y,*n).noalias() += matrix(a,*n,*n,*lda).selfadjointView() * (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; -} - -/** ZHBMV performs the matrix-vector operation - * - * y := alpha*A*x + beta*y, - * - * where alpha and beta are scalars, x and y are n element vectors and - * A is an n by n hermitian band matrix, with k super-diagonals. - */ -int EIGEN_BLAS_FUNC(hbmv)(char *uplo, int *n, int *k, RealScalar *alpha, RealScalar *a, int *lda, - RealScalar *x, int *incx, RealScalar *beta, RealScalar *y, int *incy) -{ - return 1; -} - -/** ZHPMV performs the matrix-vector operation - * - * y := alpha*A*x + beta*y, - * - * where alpha and beta are scalars, x and y are n element vectors and - * A is an n by n hermitian matrix, supplied in packed form. - */ -int EIGEN_BLAS_FUNC(hpmv)(char *uplo, int *n, RealScalar *alpha, RealScalar *ap, RealScalar *x, int *incx, RealScalar *beta, RealScalar *y, int *incy) -{ - return 1; -} - -/** ZHPR performs the hermitian rank 1 operation - * - * A := alpha*x*conjg( x' ) + A, - * - * where alpha is a real scalar, x is an n element vector and A is an - * n by n hermitian matrix, supplied in packed form. - */ -int EIGEN_BLAS_FUNC(hpr)(char *uplo, int *n, RealScalar *alpha, RealScalar *x, int *incx, RealScalar *ap) -{ - return 1; -} - -/** ZHPR2 performs the hermitian rank 2 operation - * - * A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A, - * - * where alpha is a scalar, x and y are n element vectors and A is an - * n by n hermitian matrix, supplied in packed form. - */ -int EIGEN_BLAS_FUNC(hpr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *x, int *incx, RealScalar *y, int *incy, RealScalar *ap) -{ - return 1; -} - -/** ZHER performs the hermitian rank 1 operation - * - * A := alpha*x*conjg( x' ) + A, - * - * 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 *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; -} + matrix(a,*m,*n,*lda) += alpha * vector(x_cpy,*m) * vector(y_cpy,*n).adjoint(); -/** ZHER2 performs the hermitian rank 2 operation - * - * A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A, - * - * 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 *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; } -/** ZGERU performs the rank 1 operation - * - * A := alpha*x*y' + A, - * - * where alpha is a scalar, x is an m element vector, y is an n element - * vector and A is an m by n matrix. - */ -int EIGEN_BLAS_FUNC(geru)(int *m, 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(*m<0) info = 1; - else if(*n<0) info = 2; - else if(*incx==0) info = 5; - else if(*incy==0) info = 7; - else if(*lda(px); - Scalar* y = reinterpret_cast(py); - Scalar* a = reinterpret_cast(pa); - Scalar alpha = *reinterpret_cast(palpha); - - int info = 0; - if(*m<0) info = 1; - else if(*n<0) info = 2; - else if(*incx==0) info = 5; - else if(*incy==0) info = 7; - else if(*lda