diff options
Diffstat (limited to 'eigen/blas/level2_impl.h')
-rw-r--r-- | eigen/blas/level2_impl.h | 553 |
1 files changed, 0 insertions, 553 deletions
diff --git a/eigen/blas/level2_impl.h b/eigen/blas/level2_impl.h deleted file mode 100644 index 173f40b..0000000 --- a/eigen/blas/level2_impl.h +++ /dev/null @@ -1,553 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2009-2010 Gael Guennebaud <gael.guennebaud@inria.fr> -// -// This Source Code Form is subject to the terms of the Mozilla -// Public License v. 2.0. If a copy of the MPL was not distributed -// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. - -#include "common.h" - -template<typename Index, typename Scalar, int StorageOrder, bool ConjugateLhs, bool ConjugateRhs> -struct general_matrix_vector_product_wrapper -{ - static void run(Index rows, Index cols,const Scalar *lhs, Index lhsStride, const Scalar *rhs, Index rhsIncr, Scalar* res, Index resIncr, Scalar alpha) - { - typedef internal::const_blas_data_mapper<Scalar,Index,StorageOrder> LhsMapper; - typedef internal::const_blas_data_mapper<Scalar,Index,RowMajor> RhsMapper; - - internal::general_matrix_vector_product - <Index,Scalar,LhsMapper,StorageOrder,ConjugateLhs,Scalar,RhsMapper,ConjugateRhs>::run( - rows, cols, LhsMapper(lhs, lhsStride), RhsMapper(rhs, rhsIncr), res, resIncr, alpha); - } -}; - -int EIGEN_BLAS_FUNC(gemv)(const char *opa, const int *m, const int *n, const RealScalar *palpha, - const RealScalar *pa, const int *lda, const RealScalar *pb, const int *incb, const RealScalar *pbeta, RealScalar *pc, const int *incc) -{ - typedef void (*functype)(int, int, const Scalar *, int, const Scalar *, int , Scalar *, int, Scalar); - static const functype func[4] = { - // array index: NOTR - (general_matrix_vector_product_wrapper<int,Scalar,ColMajor,false,false>::run), - // array index: TR - (general_matrix_vector_product_wrapper<int,Scalar,RowMajor,false,false>::run), - // array index: ADJ - (general_matrix_vector_product_wrapper<int,Scalar,RowMajor,Conj ,false>::run), - 0 - }; - - const Scalar* a = reinterpret_cast<const Scalar*>(pa); - const Scalar* b = reinterpret_cast<const Scalar*>(pb); - Scalar* c = reinterpret_cast<Scalar*>(pc); - Scalar alpha = *reinterpret_cast<const Scalar*>(palpha); - Scalar beta = *reinterpret_cast<const Scalar*>(pbeta); - - // check arguments - int info = 0; - if(OP(*opa)==INVALID) info = 1; - else if(*m<0) info = 2; - else if(*n<0) info = 3; - else if(*lda<std::max(1,*m)) info = 6; - else if(*incb==0) info = 8; - else if(*incc==0) info = 11; - if(info) - return xerbla_(SCALAR_SUFFIX_UP"GEMV ",&info,6); - - if(*m==0 || *n==0 || (alpha==Scalar(0) && beta==Scalar(1))) - return 0; - - int actual_m = *m; - int actual_n = *n; - int code = OP(*opa); - if(code!=NOTR) - std::swap(actual_m,actual_n); - - const Scalar* actual_b = get_compact_vector(b,actual_n,*incb); - Scalar* actual_c = get_compact_vector(c,actual_m,*incc); - - if(beta!=Scalar(1)) - { - if(beta==Scalar(0)) make_vector(actual_c, actual_m).setZero(); - else make_vector(actual_c, actual_m) *= beta; - } - - if(code>=4 || func[code]==0) - return 0; - - func[code](actual_m, actual_n, a, *lda, actual_b, 1, actual_c, 1, alpha); - - if(actual_b!=b) delete[] actual_b; - if(actual_c!=c) delete[] copy_back(actual_c,c,actual_m,*incc); - - return 1; -} - -int EIGEN_BLAS_FUNC(trsv)(const char *uplo, const char *opa, const char *diag, const int *n, const RealScalar *pa, const int *lda, RealScalar *pb, const int *incb) -{ - typedef void (*functype)(int, const Scalar *, int, Scalar *); - static const functype func[16] = { - // array index: NOTR | (UP << 2) | (NUNIT << 3) - (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0, false,ColMajor>::run), - // array index: TR | (UP << 2) | (NUNIT << 3) - (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0, false,RowMajor>::run), - // array index: ADJ | (UP << 2) | (NUNIT << 3) - (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0, Conj, RowMajor>::run), - 0, - // array index: NOTR | (LO << 2) | (NUNIT << 3) - (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0, false,ColMajor>::run), - // array index: TR | (LO << 2) | (NUNIT << 3) - (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0, false,RowMajor>::run), - // array index: ADJ | (LO << 2) | (NUNIT << 3) - (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0, Conj, RowMajor>::run), - 0, - // array index: NOTR | (UP << 2) | (UNIT << 3) - (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,ColMajor>::run), - // array index: TR | (UP << 2) | (UNIT << 3) - (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,RowMajor>::run), - // array index: ADJ | (UP << 2) | (UNIT << 3) - (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,Conj, RowMajor>::run), - 0, - // array index: NOTR | (LO << 2) | (UNIT << 3) - (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,ColMajor>::run), - // array index: TR | (LO << 2) | (UNIT << 3) - (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,RowMajor>::run), - // array index: ADJ | (LO << 2) | (UNIT << 3) - (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,Conj, RowMajor>::run), - 0 - }; - - const Scalar* a = reinterpret_cast<const Scalar*>(pa); - Scalar* b = reinterpret_cast<Scalar*>(pb); - - int info = 0; - if(UPLO(*uplo)==INVALID) info = 1; - else if(OP(*opa)==INVALID) info = 2; - else if(DIAG(*diag)==INVALID) info = 3; - else if(*n<0) info = 4; - else if(*lda<std::max(1,*n)) info = 6; - else if(*incb==0) info = 8; - if(info) - return xerbla_(SCALAR_SUFFIX_UP"TRSV ",&info,6); - - Scalar* actual_b = get_compact_vector(b,*n,*incb); - - int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3); - func[code](*n, a, *lda, actual_b); - - if(actual_b!=b) delete[] copy_back(actual_b,b,*n,*incb); - - return 0; -} - - - -int EIGEN_BLAS_FUNC(trmv)(const char *uplo, const char *opa, const char *diag, const int *n, const RealScalar *pa, const int *lda, RealScalar *pb, const int *incb) -{ - typedef void (*functype)(int, int, const Scalar *, int, const Scalar *, int, Scalar *, int, const Scalar&); - static const functype func[16] = { - // array index: NOTR | (UP << 2) | (NUNIT << 3) - (internal::triangular_matrix_vector_product<int,Upper|0, Scalar,false,Scalar,false,ColMajor>::run), - // array index: TR | (UP << 2) | (NUNIT << 3) - (internal::triangular_matrix_vector_product<int,Lower|0, Scalar,false,Scalar,false,RowMajor>::run), - // array index: ADJ | (UP << 2) | (NUNIT << 3) - (internal::triangular_matrix_vector_product<int,Lower|0, Scalar,Conj, Scalar,false,RowMajor>::run), - 0, - // array index: NOTR | (LO << 2) | (NUNIT << 3) - (internal::triangular_matrix_vector_product<int,Lower|0, Scalar,false,Scalar,false,ColMajor>::run), - // array index: TR | (LO << 2) | (NUNIT << 3) - (internal::triangular_matrix_vector_product<int,Upper|0, Scalar,false,Scalar,false,RowMajor>::run), - // array index: ADJ | (LO << 2) | (NUNIT << 3) - (internal::triangular_matrix_vector_product<int,Upper|0, Scalar,Conj, Scalar,false,RowMajor>::run), - 0, - // array index: NOTR | (UP << 2) | (UNIT << 3) - (internal::triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run), - // array index: TR | (UP << 2) | (UNIT << 3) - (internal::triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run), - // array index: ADJ | (UP << 2) | (UNIT << 3) - (internal::triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run), - 0, - // array index: NOTR | (LO << 2) | (UNIT << 3) - (internal::triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run), - // array index: TR | (LO << 2) | (UNIT << 3) - (internal::triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run), - // array index: ADJ | (LO << 2) | (UNIT << 3) - (internal::triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run), - 0 - }; - - const Scalar* a = reinterpret_cast<const Scalar*>(pa); - Scalar* b = reinterpret_cast<Scalar*>(pb); - - int info = 0; - if(UPLO(*uplo)==INVALID) info = 1; - else if(OP(*opa)==INVALID) info = 2; - else if(DIAG(*diag)==INVALID) info = 3; - else if(*n<0) info = 4; - else if(*lda<std::max(1,*n)) info = 6; - else if(*incb==0) info = 8; - if(info) - return xerbla_(SCALAR_SUFFIX_UP"TRMV ",&info,6); - - if(*n==0) - return 1; - - Scalar* actual_b = get_compact_vector(b,*n,*incb); - Matrix<Scalar,Dynamic,1> res(*n); - res.setZero(); - - int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3); - if(code>=16 || func[code]==0) - return 0; - - func[code](*n, *n, a, *lda, actual_b, 1, res.data(), 1, Scalar(1)); - - copy_back(res.data(),b,*n,*incb); - if(actual_b!=b) delete[] actual_b; - - return 1; -} - -/** GBMV performs one of the matrix-vector operations - * - * y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, - * - * 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 *palpha, RealScalar *pa, int *lda, - RealScalar *px, int *incx, RealScalar *pbeta, RealScalar *py, int *incy) -{ - const Scalar* a = reinterpret_cast<const Scalar*>(pa); - const Scalar* x = reinterpret_cast<const Scalar*>(px); - Scalar* y = reinterpret_cast<Scalar*>(py); - Scalar alpha = *reinterpret_cast<const Scalar*>(palpha); - Scalar beta = *reinterpret_cast<const Scalar*>(pbeta); - int coeff_rows = *kl+*ku+1; - - int info = 0; - if(OP(*trans)==INVALID) info = 1; - else if(*m<0) info = 2; - else if(*n<0) info = 3; - else if(*kl<0) info = 4; - else if(*ku<0) info = 5; - else if(*lda<coeff_rows) info = 8; - else if(*incx==0) info = 10; - else if(*incy==0) info = 13; - if(info) - return xerbla_(SCALAR_SUFFIX_UP"GBMV ",&info,6); - - if(*m==0 || *n==0 || (alpha==Scalar(0) && beta==Scalar(1))) - return 0; - - int actual_m = *m; - int actual_n = *n; - if(OP(*trans)!=NOTR) - std::swap(actual_m,actual_n); - - const Scalar* actual_x = get_compact_vector(x,actual_n,*incx); - Scalar* actual_y = get_compact_vector(y,actual_m,*incy); - - if(beta!=Scalar(1)) - { - if(beta==Scalar(0)) make_vector(actual_y, actual_m).setZero(); - else make_vector(actual_y, actual_m) *= beta; - } - - ConstMatrixType mat_coeffs(a,coeff_rows,*n,*lda); - - int nb = std::min(*n,(*m)+(*ku)); - for(int j=0; j<nb; ++j) - { - int start = std::max(0,j - *ku); - int end = std::min((*m)-1,j + *kl); - int len = end - start + 1; - int offset = (*ku) - j + start; - if(OP(*trans)==NOTR) - make_vector(actual_y+start,len) += (alpha*actual_x[j]) * mat_coeffs.col(j).segment(offset,len); - else if(OP(*trans)==TR) - actual_y[j] += alpha * ( mat_coeffs.col(j).segment(offset,len).transpose() * make_vector(actual_x+start,len) ).value(); - else - actual_y[j] += alpha * ( mat_coeffs.col(j).segment(offset,len).adjoint() * make_vector(actual_x+start,len) ).value(); - } - - if(actual_x!=x) delete[] actual_x; - if(actual_y!=y) delete[] copy_back(actual_y,y,actual_m,*incy); - - return 0; -} - -#if 0 -/** TBMV performs one of the matrix-vector operations - * - * x := A*x, or x := A'*x, - * - * 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 *opa, char *diag, int *n, int *k, RealScalar *pa, int *lda, RealScalar *px, int *incx) -{ - Scalar* a = reinterpret_cast<Scalar*>(pa); - Scalar* x = reinterpret_cast<Scalar*>(px); - int coeff_rows = *k + 1; - - int info = 0; - if(UPLO(*uplo)==INVALID) info = 1; - else if(OP(*opa)==INVALID) info = 2; - else if(DIAG(*diag)==INVALID) info = 3; - else if(*n<0) info = 4; - else if(*k<0) info = 5; - else if(*lda<coeff_rows) info = 7; - else if(*incx==0) info = 9; - if(info) - return xerbla_(SCALAR_SUFFIX_UP"TBMV ",&info,6); - - if(*n==0) - return 0; - - int actual_n = *n; - - Scalar* actual_x = get_compact_vector(x,actual_n,*incx); - - MatrixType mat_coeffs(a,coeff_rows,*n,*lda); - - int ku = UPLO(*uplo)==UPPER ? *k : 0; - int kl = UPLO(*uplo)==LOWER ? *k : 0; - - for(int j=0; j<*n; ++j) - { - int start = std::max(0,j - ku); - int end = std::min((*m)-1,j + kl); - int len = end - start + 1; - int offset = (ku) - j + start; - - if(OP(*trans)==NOTR) - make_vector(actual_y+start,len) += (alpha*actual_x[j]) * mat_coeffs.col(j).segment(offset,len); - else if(OP(*trans)==TR) - actual_y[j] += alpha * ( mat_coeffs.col(j).segment(offset,len).transpose() * make_vector(actual_x+start,len) ).value(); - else - actual_y[j] += alpha * ( mat_coeffs.col(j).segment(offset,len).adjoint() * make_vector(actual_x+start,len) ).value(); - } - - if(actual_x!=x) delete[] actual_x; - if(actual_y!=y) delete[] copy_back(actual_y,y,actual_m,*incy); - - return 0; -} -#endif - -/** DTBSV solves one of the systems of equations - * - * A*x = b, or A'*x = b, - * - * where b and x are n element vectors and A is an n by n unit, or - * non-unit, upper or lower triangular band matrix, with ( k + 1 ) - * diagonals. - * - * 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 *op, char *diag, int *n, int *k, RealScalar *pa, int *lda, RealScalar *px, int *incx) -{ - typedef void (*functype)(int, int, const Scalar *, int, Scalar *); - static const functype func[16] = { - // array index: NOTR | (UP << 2) | (NUNIT << 3) - (internal::band_solve_triangular_selector<int,Upper|0, Scalar,false,Scalar,ColMajor>::run), - // array index: TR | (UP << 2) | (NUNIT << 3) - (internal::band_solve_triangular_selector<int,Lower|0, Scalar,false,Scalar,RowMajor>::run), - // array index: ADJ | (UP << 2) | (NUNIT << 3) - (internal::band_solve_triangular_selector<int,Lower|0, Scalar,Conj, Scalar,RowMajor>::run), - 0, - // array index: NOTR | (LO << 2) | (NUNIT << 3) - (internal::band_solve_triangular_selector<int,Lower|0, Scalar,false,Scalar,ColMajor>::run), - // array index: TR | (LO << 2) | (NUNIT << 3) - (internal::band_solve_triangular_selector<int,Upper|0, Scalar,false,Scalar,RowMajor>::run), - // array index: ADJ | (LO << 2) | (NUNIT << 3) - (internal::band_solve_triangular_selector<int,Upper|0, Scalar,Conj, Scalar,RowMajor>::run), - 0, - // array index: NOTR | (UP << 2) | (UNIT << 3) - (internal::band_solve_triangular_selector<int,Upper|UnitDiag,Scalar,false,Scalar,ColMajor>::run), - // array index: TR | (UP << 2) | (UNIT << 3) - (internal::band_solve_triangular_selector<int,Lower|UnitDiag,Scalar,false,Scalar,RowMajor>::run), - // array index: ADJ | (UP << 2) | (UNIT << 3) - (internal::band_solve_triangular_selector<int,Lower|UnitDiag,Scalar,Conj, Scalar,RowMajor>::run), - 0, - // array index: NOTR | (LO << 2) | (UNIT << 3) - (internal::band_solve_triangular_selector<int,Lower|UnitDiag,Scalar,false,Scalar,ColMajor>::run), - // array index: TR | (LO << 2) | (UNIT << 3) - (internal::band_solve_triangular_selector<int,Upper|UnitDiag,Scalar,false,Scalar,RowMajor>::run), - // array index: ADJ | (LO << 2) | (UNIT << 3) - (internal::band_solve_triangular_selector<int,Upper|UnitDiag,Scalar,Conj, Scalar,RowMajor>::run), - 0, - }; - - Scalar* a = reinterpret_cast<Scalar*>(pa); - Scalar* x = reinterpret_cast<Scalar*>(px); - int coeff_rows = *k+1; - - int info = 0; - if(UPLO(*uplo)==INVALID) info = 1; - else if(OP(*op)==INVALID) info = 2; - else if(DIAG(*diag)==INVALID) info = 3; - else if(*n<0) info = 4; - else if(*k<0) info = 5; - else if(*lda<coeff_rows) info = 7; - else if(*incx==0) info = 9; - if(info) - return xerbla_(SCALAR_SUFFIX_UP"TBSV ",&info,6); - - if(*n==0 || (*k==0 && DIAG(*diag)==UNIT)) - return 0; - - int actual_n = *n; - - Scalar* actual_x = get_compact_vector(x,actual_n,*incx); - - int code = OP(*op) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3); - if(code>=16 || func[code]==0) - return 0; - - func[code](*n, *k, a, *lda, actual_x); - - if(actual_x!=x) delete[] copy_back(actual_x,x,actual_n,*incx); - - return 0; -} - -/** DTPMV performs one of the matrix-vector operations - * - * x := A*x, or x := A'*x, - * - * 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 *opa, char *diag, int *n, RealScalar *pap, RealScalar *px, int *incx) -{ - typedef void (*functype)(int, const Scalar*, const Scalar*, Scalar*, Scalar); - static const functype func[16] = { - // array index: NOTR | (UP << 2) | (NUNIT << 3) - (internal::packed_triangular_matrix_vector_product<int,Upper|0, Scalar,false,Scalar,false,ColMajor>::run), - // array index: TR | (UP << 2) | (NUNIT << 3) - (internal::packed_triangular_matrix_vector_product<int,Lower|0, Scalar,false,Scalar,false,RowMajor>::run), - // array index: ADJ | (UP << 2) | (NUNIT << 3) - (internal::packed_triangular_matrix_vector_product<int,Lower|0, Scalar,Conj, Scalar,false,RowMajor>::run), - 0, - // array index: NOTR | (LO << 2) | (NUNIT << 3) - (internal::packed_triangular_matrix_vector_product<int,Lower|0, Scalar,false,Scalar,false,ColMajor>::run), - // array index: TR | (LO << 2) | (NUNIT << 3) - (internal::packed_triangular_matrix_vector_product<int,Upper|0, Scalar,false,Scalar,false,RowMajor>::run), - // array index: ADJ | (LO << 2) | (NUNIT << 3) - (internal::packed_triangular_matrix_vector_product<int,Upper|0, Scalar,Conj, Scalar,false,RowMajor>::run), - 0, - // array index: NOTR | (UP << 2) | (UNIT << 3) - (internal::packed_triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run), - // array index: TR | (UP << 2) | (UNIT << 3) - (internal::packed_triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run), - // array index: ADJ | (UP << 2) | (UNIT << 3) - (internal::packed_triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run), - 0, - // array index: NOTR | (LO << 2) | (UNIT << 3) - (internal::packed_triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run), - // array index: TR | (LO << 2) | (UNIT << 3) - (internal::packed_triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run), - // array index: ADJ | (LO << 2) | (UNIT << 3) - (internal::packed_triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run), - 0 - }; - - Scalar* ap = reinterpret_cast<Scalar*>(pap); - Scalar* x = reinterpret_cast<Scalar*>(px); - - int info = 0; - if(UPLO(*uplo)==INVALID) info = 1; - else if(OP(*opa)==INVALID) info = 2; - else if(DIAG(*diag)==INVALID) info = 3; - else if(*n<0) info = 4; - else if(*incx==0) info = 7; - if(info) - return xerbla_(SCALAR_SUFFIX_UP"TPMV ",&info,6); - - if(*n==0) - return 1; - - Scalar* actual_x = get_compact_vector(x,*n,*incx); - Matrix<Scalar,Dynamic,1> res(*n); - res.setZero(); - - int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3); - if(code>=16 || func[code]==0) - return 0; - - func[code](*n, ap, actual_x, res.data(), Scalar(1)); - - copy_back(res.data(),x,*n,*incx); - if(actual_x!=x) delete[] actual_x; - - return 1; -} - -/** DTPSV solves one of the systems of equations - * - * A*x = b, or A'*x = b, - * - * where b and x are n element vectors and A is an n by n unit, or - * non-unit, upper or lower triangular matrix, supplied in packed form. - * - * 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 *opa, char *diag, int *n, RealScalar *pap, RealScalar *px, int *incx) -{ - typedef void (*functype)(int, const Scalar*, Scalar*); - static const functype func[16] = { - // array index: NOTR | (UP << 2) | (NUNIT << 3) - (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0, false,ColMajor>::run), - // array index: TR | (UP << 2) | (NUNIT << 3) - (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0, false,RowMajor>::run), - // array index: ADJ | (UP << 2) | (NUNIT << 3) - (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0, Conj, RowMajor>::run), - 0, - // array index: NOTR | (LO << 2) | (NUNIT << 3) - (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0, false,ColMajor>::run), - // array index: TR | (LO << 2) | (NUNIT << 3) - (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0, false,RowMajor>::run), - // array index: ADJ | (LO << 2) | (NUNIT << 3) - (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0, Conj, RowMajor>::run), - 0, - // array index: NOTR | (UP << 2) | (UNIT << 3) - (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,ColMajor>::run), - // array index: TR | (UP << 2) | (UNIT << 3) - (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,RowMajor>::run), - // array index: ADJ | (UP << 2) | (UNIT << 3) - (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,Conj, RowMajor>::run), - 0, - // array index: NOTR | (LO << 2) | (UNIT << 3) - (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,ColMajor>::run), - // array index: TR | (LO << 2) | (UNIT << 3) - (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,RowMajor>::run), - // array index: ADJ | (LO << 2) | (UNIT << 3) - (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,Conj, RowMajor>::run), - 0 - }; - - Scalar* ap = reinterpret_cast<Scalar*>(pap); - Scalar* x = reinterpret_cast<Scalar*>(px); - - int info = 0; - if(UPLO(*uplo)==INVALID) info = 1; - else if(OP(*opa)==INVALID) info = 2; - else if(DIAG(*diag)==INVALID) info = 3; - else if(*n<0) info = 4; - else if(*incx==0) info = 7; - if(info) - return xerbla_(SCALAR_SUFFIX_UP"TPSV ",&info,6); - - Scalar* actual_x = get_compact_vector(x,*n,*incx); - - int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3); - func[code](*n, ap, actual_x); - - if(actual_x!=x) delete[] copy_back(actual_x,x,*n,*incx); - - return 1; -} |