| // 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); |
| } |
| }; |
| |
| 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); |
| |
| if (*m == 0 || *n == 0 || (alpha == Scalar(0) && beta == Scalar(1))) return; |
| |
| 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; |
| |
| 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); |
| } |
| |
| 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); |
| |
| 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); |
| } |
| |
| 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); |
| |
| if (*n == 0) return; |
| |
| 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; |
| |
| 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; |
| } |
| |
| /** 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. |
| */ |
| 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); |
| |
| if (*m == 0 || *n == 0 || (alpha == Scalar(0) && beta == Scalar(1))) return; |
| |
| 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); |
| } |
| |
| #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. |
| */ |
| 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; |
| |
| 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); |
| } |
| #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. |
| */ |
| 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); |
| |
| if (*n == 0 || (*k == 0 && DIAG(*diag) == UNIT)) return; |
| |
| 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; |
| |
| func[code](*n, *k, a, *lda, actual_x); |
| |
| if (actual_x != x) delete[] copy_back(actual_x, x, actual_n, *incx); |
| } |
| |
| /** 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. |
| */ |
| 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); |
| |
| if (*n == 0) return; |
| |
| 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; |
| |
| func[code](*n, ap, actual_x, res.data(), Scalar(1)); |
| |
| copy_back(res.data(), x, *n, *incx); |
| if (actual_x != x) delete[] actual_x; |
| } |
| |
| /** 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. |
| */ |
| 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); |
| |
| 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); |
| } |
