| // 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/. | |
| // SPDX-License-Identifier: 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 Eigen::internal::const_blas_data_mapper<Scalar, Index, StorageOrder> LhsMapper; | |
| typedef Eigen::internal::const_blas_data_mapper<Scalar, Index, Eigen::RowMajor> RhsMapper; | |
| Eigen::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 EIGEN_BLAS_INT *m, const EIGEN_BLAS_INT *n, const RealScalar *palpha, const RealScalar *pa, | |
| const EIGEN_BLAS_INT *lda, const RealScalar *pb, const EIGEN_BLAS_INT *incb, const RealScalar *pbeta, RealScalar *pc, | |
| const EIGEN_BLAS_INT *incc) { | |
| typedef void (*functype)(EIGEN_BLAS_INT, EIGEN_BLAS_INT, const Scalar *, EIGEN_BLAS_INT, const Scalar *, | |
| EIGEN_BLAS_INT, Scalar *, EIGEN_BLAS_INT, Scalar); | |
| static const functype func[4] = { | |
| // array index: NOTR | |
| (general_matrix_vector_product_wrapper<EIGEN_BLAS_INT, Scalar, Eigen::ColMajor, false, false>::run), | |
| // array index: TR | |
| (general_matrix_vector_product_wrapper<EIGEN_BLAS_INT, Scalar, Eigen::RowMajor, false, false>::run), | |
| // array index: ADJ | |
| (general_matrix_vector_product_wrapper<EIGEN_BLAS_INT, Scalar, Eigen::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 | |
| EIGEN_BLAS_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<EIGEN_BLAS_INT>(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, kBlasNameLength); | |
| if (*m == 0 || *n == 0 || (alpha == Scalar(0) && beta == Scalar(1))) return; | |
| EIGEN_BLAS_INT actual_m = *m; | |
| EIGEN_BLAS_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 EIGEN_BLAS_INT *n, const RealScalar *pa, | |
| const EIGEN_BLAS_INT *lda, RealScalar *pb, const EIGEN_BLAS_INT *incb) { | |
| typedef void (*functype)(EIGEN_BLAS_INT, const Scalar *, EIGEN_BLAS_INT, Scalar *); | |
| using Eigen::ColMajor; | |
| using Eigen::Lower; | |
| using Eigen::OnTheLeft; | |
| using Eigen::RowMajor; | |
| using Eigen::UnitDiag; | |
| using Eigen::Upper; | |
| static const functype func[16] = {// array index: NOTR | (UP << 2) | (NUNIT << 3) | |
| (Eigen::internal::triangular_solve_vector<Scalar, Scalar, EIGEN_BLAS_INT, OnTheLeft, | |
| Upper | 0, false, ColMajor>::run), | |
| // array index: TR | (UP << 2) | (NUNIT << 3) | |
| (Eigen::internal::triangular_solve_vector<Scalar, Scalar, EIGEN_BLAS_INT, OnTheLeft, | |
| Lower | 0, false, RowMajor>::run), | |
| // array index: ADJ | (UP << 2) | (NUNIT << 3) | |
| (Eigen::internal::triangular_solve_vector<Scalar, Scalar, EIGEN_BLAS_INT, OnTheLeft, | |
| Lower | 0, Conj, RowMajor>::run), | |
| 0, | |
| // array index: NOTR | (LO << 2) | (NUNIT << 3) | |
| (Eigen::internal::triangular_solve_vector<Scalar, Scalar, EIGEN_BLAS_INT, OnTheLeft, | |
| Lower | 0, false, ColMajor>::run), | |
| // array index: TR | (LO << 2) | (NUNIT << 3) | |
| (Eigen::internal::triangular_solve_vector<Scalar, Scalar, EIGEN_BLAS_INT, OnTheLeft, | |
| Upper | 0, false, RowMajor>::run), | |
| // array index: ADJ | (LO << 2) | (NUNIT << 3) | |
| (Eigen::internal::triangular_solve_vector<Scalar, Scalar, EIGEN_BLAS_INT, OnTheLeft, | |
| Upper | 0, Conj, RowMajor>::run), | |
| 0, | |
| // array index: NOTR | (UP << 2) | (UNIT << 3) | |
| (Eigen::internal::triangular_solve_vector<Scalar, Scalar, EIGEN_BLAS_INT, OnTheLeft, | |
| Upper | UnitDiag, false, ColMajor>::run), | |
| // array index: TR | (UP << 2) | (UNIT << 3) | |
| (Eigen::internal::triangular_solve_vector<Scalar, Scalar, EIGEN_BLAS_INT, OnTheLeft, | |
| Lower | UnitDiag, false, RowMajor>::run), | |
| // array index: ADJ | (UP << 2) | (UNIT << 3) | |
| (Eigen::internal::triangular_solve_vector<Scalar, Scalar, EIGEN_BLAS_INT, OnTheLeft, | |
| Lower | UnitDiag, Conj, RowMajor>::run), | |
| 0, | |
| // array index: NOTR | (LO << 2) | (UNIT << 3) | |
| (Eigen::internal::triangular_solve_vector<Scalar, Scalar, EIGEN_BLAS_INT, OnTheLeft, | |
| Lower | UnitDiag, false, ColMajor>::run), | |
| // array index: TR | (LO << 2) | (UNIT << 3) | |
| (Eigen::internal::triangular_solve_vector<Scalar, Scalar, EIGEN_BLAS_INT, OnTheLeft, | |
| Upper | UnitDiag, false, RowMajor>::run), | |
| // array index: ADJ | (LO << 2) | (UNIT << 3) | |
| (Eigen::internal::triangular_solve_vector<Scalar, Scalar, EIGEN_BLAS_INT, OnTheLeft, | |
| Upper | UnitDiag, Conj, RowMajor>::run), | |
| 0}; | |
| const Scalar *a = reinterpret_cast<const Scalar *>(pa); | |
| Scalar *b = reinterpret_cast<Scalar *>(pb); | |
| EIGEN_BLAS_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<EIGEN_BLAS_INT>(1, *n)) | |
| info = 6; | |
| else if (*incb == 0) | |
| info = 8; | |
| if (info) return xerbla_(SCALAR_SUFFIX_UP "TRSV ", &info, kBlasNameLength); | |
| 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 EIGEN_BLAS_INT *n, const RealScalar *pa, | |
| const EIGEN_BLAS_INT *lda, RealScalar *pb, const EIGEN_BLAS_INT *incb) { | |
| typedef void (*functype)(EIGEN_BLAS_INT, EIGEN_BLAS_INT, const Scalar *, EIGEN_BLAS_INT, const Scalar *, | |
| EIGEN_BLAS_INT, Scalar *, EIGEN_BLAS_INT, const Scalar &); | |
| using Eigen::ColMajor; | |
| using Eigen::Lower; | |
| using Eigen::OnTheLeft; | |
| using Eigen::RowMajor; | |
| using Eigen::UnitDiag; | |
| using Eigen::Upper; | |
| static const functype func[16] = { | |
| // array index: NOTR | (UP << 2) | (NUNIT << 3) | |
| (Eigen::internal::triangular_matrix_vector_product<EIGEN_BLAS_INT, Upper | 0, Scalar, false, Scalar, false, | |
| ColMajor>::run), | |
| // array index: TR | (UP << 2) | (NUNIT << 3) | |
| (Eigen::internal::triangular_matrix_vector_product<EIGEN_BLAS_INT, Lower | 0, Scalar, false, Scalar, false, | |
| RowMajor>::run), | |
| // array index: ADJ | (UP << 2) | (NUNIT << 3) | |
| (Eigen::internal::triangular_matrix_vector_product<EIGEN_BLAS_INT, Lower | 0, Scalar, Conj, Scalar, false, | |
| RowMajor>::run), | |
| 0, | |
| // array index: NOTR | (LO << 2) | (NUNIT << 3) | |
| (Eigen::internal::triangular_matrix_vector_product<EIGEN_BLAS_INT, Lower | 0, Scalar, false, Scalar, false, | |
| ColMajor>::run), | |
| // array index: TR | (LO << 2) | (NUNIT << 3) | |
| (Eigen::internal::triangular_matrix_vector_product<EIGEN_BLAS_INT, Upper | 0, Scalar, false, Scalar, false, | |
| RowMajor>::run), | |
| // array index: ADJ | (LO << 2) | (NUNIT << 3) | |
| (Eigen::internal::triangular_matrix_vector_product<EIGEN_BLAS_INT, Upper | 0, Scalar, Conj, Scalar, false, | |
| RowMajor>::run), | |
| 0, | |
| // array index: NOTR | (UP << 2) | (UNIT << 3) | |
| (Eigen::internal::triangular_matrix_vector_product<EIGEN_BLAS_INT, Upper | UnitDiag, Scalar, false, Scalar, false, | |
| ColMajor>::run), | |
| // array index: TR | (UP << 2) | (UNIT << 3) | |
| (Eigen::internal::triangular_matrix_vector_product<EIGEN_BLAS_INT, Lower | UnitDiag, Scalar, false, Scalar, false, | |
| RowMajor>::run), | |
| // array index: ADJ | (UP << 2) | (UNIT << 3) | |
| (Eigen::internal::triangular_matrix_vector_product<EIGEN_BLAS_INT, Lower | UnitDiag, Scalar, Conj, Scalar, false, | |
| RowMajor>::run), | |
| 0, | |
| // array index: NOTR | (LO << 2) | (UNIT << 3) | |
| (Eigen::internal::triangular_matrix_vector_product<EIGEN_BLAS_INT, Lower | UnitDiag, Scalar, false, Scalar, false, | |
| ColMajor>::run), | |
| // array index: TR | (LO << 2) | (UNIT << 3) | |
| (Eigen::internal::triangular_matrix_vector_product<EIGEN_BLAS_INT, Upper | UnitDiag, Scalar, false, Scalar, false, | |
| RowMajor>::run), | |
| // array index: ADJ | (LO << 2) | (UNIT << 3) | |
| (Eigen::internal::triangular_matrix_vector_product<EIGEN_BLAS_INT, Upper | UnitDiag, Scalar, Conj, Scalar, false, | |
| RowMajor>::run), | |
| 0}; | |
| const Scalar *a = reinterpret_cast<const Scalar *>(pa); | |
| Scalar *b = reinterpret_cast<Scalar *>(pb); | |
| EIGEN_BLAS_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<EIGEN_BLAS_INT>(1, *n)) | |
| info = 6; | |
| else if (*incb == 0) | |
| info = 8; | |
| if (info) return xerbla_(SCALAR_SUFFIX_UP "TRMV ", &info, kBlasNameLength); | |
| if (*n == 0) return; | |
| Scalar *actual_b = get_compact_vector(b, *n, *incb); | |
| Eigen::Matrix<Scalar, Eigen::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, EIGEN_BLAS_INT *m, EIGEN_BLAS_INT *n, EIGEN_BLAS_INT *kl, EIGEN_BLAS_INT *ku, RealScalar *palpha, | |
| RealScalar *pa, EIGEN_BLAS_INT *lda, RealScalar *px, EIGEN_BLAS_INT *incx, RealScalar *pbeta, RealScalar *py, | |
| EIGEN_BLAS_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); | |
| EIGEN_BLAS_INT coeff_rows = *kl + *ku + 1; | |
| EIGEN_BLAS_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, kBlasNameLength); | |
| if (*m == 0 || *n == 0 || (alpha == Scalar(0) && beta == Scalar(1))) return; | |
| EIGEN_BLAS_INT actual_m = *m; | |
| EIGEN_BLAS_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); | |
| EIGEN_BLAS_INT nb = alpha == Scalar(0) ? 0 : std::min(*n, (*m) + (*ku)); | |
| for (EIGEN_BLAS_INT j = 0; j < nb; ++j) { | |
| EIGEN_BLAS_INT start = std::max<EIGEN_BLAS_INT>(0, j - *ku); | |
| EIGEN_BLAS_INT end = std::min((*m) - 1, j + *kl); | |
| EIGEN_BLAS_INT len = end - start + 1; | |
| EIGEN_BLAS_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); | |
| } | |
| /** TBMV performs one of the matrix-vector operations | |
| * | |
| * x := A*x, or x := A'*x, or x := conjg(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. | |
| * | |
| * Band storage: upper triangle stores A[i,j] at a[(k+i-j) + j*lda], | |
| * lower triangle stores A[i,j] at a[(i-j) + j*lda]. | |
| */ | |
| EIGEN_BLAS_FUNC(tbmv) | |
| (char *uplo, char *opa, char *diag, EIGEN_BLAS_INT *n, EIGEN_BLAS_INT *k, RealScalar *pa, EIGEN_BLAS_INT *lda, | |
| RealScalar *px, EIGEN_BLAS_INT *incx) { | |
| Scalar *a = reinterpret_cast<Scalar *>(pa); | |
| Scalar *x = reinterpret_cast<Scalar *>(px); | |
| EIGEN_BLAS_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 < *k + 1) | |
| info = 7; | |
| else if (*incx == 0) | |
| info = 9; | |
| if (info) return xerbla_(SCALAR_SUFFIX_UP "TBMV ", &info, kBlasNameLength); | |
| if (*n == 0) return; | |
| Scalar *actual_x = get_compact_vector(x, *n, *incx); | |
| bool upper = (UPLO(*uplo) == UP); | |
| int op = OP(*opa); | |
| bool unit = (DIAG(*diag) == UNIT); | |
| if (*k >= 8) { | |
| // Vectorized path: use Eigen Map segments for the inner band operations. | |
| ConstMatrixType band(a, *k + 1, *n, *lda); | |
| if (op == NOTR) { | |
| if (upper) { | |
| for (EIGEN_BLAS_INT j = 0; j < *n; ++j) { | |
| if (actual_x[j] != Scalar(0)) { | |
| EIGEN_BLAS_INT start = std::max<EIGEN_BLAS_INT>(0, j - *k); | |
| EIGEN_BLAS_INT len = j - start; | |
| EIGEN_BLAS_INT offset = *k - (j - start); | |
| Scalar temp = actual_x[j]; | |
| if (len > 0) make_vector(actual_x + start, len) += temp * band.col(j).segment(offset, len); | |
| if (!unit) actual_x[j] = temp * band(*k, j); | |
| } | |
| } | |
| } else { | |
| for (EIGEN_BLAS_INT j = *n - 1; j >= 0; --j) { | |
| if (actual_x[j] != Scalar(0)) { | |
| EIGEN_BLAS_INT len = std::min(*n - 1, j + *k) - j; | |
| Scalar temp = actual_x[j]; | |
| if (len > 0) make_vector(actual_x + j + 1, len) += temp * band.col(j).segment(1, len); | |
| if (!unit) actual_x[j] = temp * band(0, j); | |
| } | |
| } | |
| } | |
| } else if (op == TR) { | |
| if (upper) { | |
| for (EIGEN_BLAS_INT j = *n - 1; j >= 0; --j) { | |
| EIGEN_BLAS_INT start = std::max<EIGEN_BLAS_INT>(0, j - *k); | |
| EIGEN_BLAS_INT len = j - start; | |
| EIGEN_BLAS_INT offset = *k - (j - start); | |
| Scalar temp = actual_x[j]; | |
| if (!unit) temp *= band(*k, j); | |
| if (len > 0) | |
| temp += (band.col(j).segment(offset, len).cwiseProduct(make_vector(actual_x + start, len))).sum(); | |
| actual_x[j] = temp; | |
| } | |
| } else { | |
| for (EIGEN_BLAS_INT j = 0; j < *n; ++j) { | |
| EIGEN_BLAS_INT len = std::min(*n - 1, j + *k) - j; | |
| Scalar temp = actual_x[j]; | |
| if (!unit) temp *= band(0, j); | |
| if (len > 0) temp += (band.col(j).segment(1, len).cwiseProduct(make_vector(actual_x + j + 1, len))).sum(); | |
| actual_x[j] = temp; | |
| } | |
| } | |
| } else { | |
| // Conjugate transpose: .dot() computes conj(lhs) . rhs. | |
| if (upper) { | |
| for (EIGEN_BLAS_INT j = *n - 1; j >= 0; --j) { | |
| EIGEN_BLAS_INT start = std::max<EIGEN_BLAS_INT>(0, j - *k); | |
| EIGEN_BLAS_INT len = j - start; | |
| EIGEN_BLAS_INT offset = *k - (j - start); | |
| Scalar temp = actual_x[j]; | |
| if (!unit) temp *= Eigen::numext::conj(band(*k, j)); | |
| if (len > 0) temp += band.col(j).segment(offset, len).dot(make_vector(actual_x + start, len)); | |
| actual_x[j] = temp; | |
| } | |
| } else { | |
| for (EIGEN_BLAS_INT j = 0; j < *n; ++j) { | |
| EIGEN_BLAS_INT len = std::min(*n - 1, j + *k) - j; | |
| Scalar temp = actual_x[j]; | |
| if (!unit) temp *= Eigen::numext::conj(band(0, j)); | |
| if (len > 0) temp += band.col(j).segment(1, len).dot(make_vector(actual_x + j + 1, len)); | |
| actual_x[j] = temp; | |
| } | |
| } | |
| } | |
| } else { | |
| // Scalar path: for narrow bandwidth, avoid Map overhead. | |
| if (op == NOTR) { | |
| if (upper) { | |
| for (EIGEN_BLAS_INT j = 0; j < *n; ++j) { | |
| if (actual_x[j] != Scalar(0)) { | |
| Scalar temp = actual_x[j]; | |
| for (EIGEN_BLAS_INT i = std::max<EIGEN_BLAS_INT>(0, j - *k); i < j; ++i) | |
| actual_x[i] += temp * a[(*k + i - j) + j * *lda]; | |
| if (!unit) actual_x[j] = temp * a[*k + j * *lda]; | |
| } | |
| } | |
| } else { | |
| for (EIGEN_BLAS_INT j = *n - 1; j >= 0; --j) { | |
| if (actual_x[j] != Scalar(0)) { | |
| Scalar temp = actual_x[j]; | |
| for (EIGEN_BLAS_INT i = j + 1; i <= std::min(*n - 1, j + *k); ++i) | |
| actual_x[i] += temp * a[(i - j) + j * *lda]; | |
| if (!unit) actual_x[j] = temp * a[j * *lda]; | |
| } | |
| } | |
| } | |
| } else { | |
| // Transpose or conjugate transpose. | |
| auto maybe_conj = [op](Scalar val) -> Scalar { return op == ADJ ? Eigen::numext::conj(val) : val; }; | |
| if (upper) { | |
| for (EIGEN_BLAS_INT j = *n - 1; j >= 0; --j) { | |
| Scalar temp = actual_x[j]; | |
| if (!unit) temp *= maybe_conj(a[*k + j * *lda]); | |
| for (EIGEN_BLAS_INT i = std::max<EIGEN_BLAS_INT>(0, j - *k); i < j; ++i) | |
| temp += maybe_conj(a[(*k + i - j) + j * *lda]) * actual_x[i]; | |
| actual_x[j] = temp; | |
| } | |
| } else { | |
| for (EIGEN_BLAS_INT j = 0; j < *n; ++j) { | |
| Scalar temp = actual_x[j]; | |
| if (!unit) temp *= maybe_conj(a[j * *lda]); | |
| for (EIGEN_BLAS_INT i = j + 1; i <= std::min(*n - 1, j + *k); ++i) | |
| temp += maybe_conj(a[(i - j) + j * *lda]) * actual_x[i]; | |
| actual_x[j] = temp; | |
| } | |
| } | |
| } | |
| } | |
| if (actual_x != x) delete[] copy_back(actual_x, x, *n, *incx); | |
| } | |
| /** 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, EIGEN_BLAS_INT *n, EIGEN_BLAS_INT *k, RealScalar *pa, EIGEN_BLAS_INT *lda, | |
| RealScalar *px, EIGEN_BLAS_INT *incx) { | |
| typedef void (*functype)(EIGEN_BLAS_INT, EIGEN_BLAS_INT, const Scalar *, EIGEN_BLAS_INT, Scalar *); | |
| using Eigen::ColMajor; | |
| using Eigen::Lower; | |
| using Eigen::OnTheLeft; | |
| using Eigen::RowMajor; | |
| using Eigen::UnitDiag; | |
| using Eigen::Upper; | |
| static const functype func[16] = { | |
| // array index: NOTR | (UP << 2) | (NUNIT << 3) | |
| (Eigen::internal::band_solve_triangular_selector<EIGEN_BLAS_INT, Upper | 0, Scalar, false, Scalar, | |
| ColMajor>::run), | |
| // array index: TR | (UP << 2) | (NUNIT << 3) | |
| (Eigen::internal::band_solve_triangular_selector<EIGEN_BLAS_INT, Lower | 0, Scalar, false, Scalar, | |
| RowMajor>::run), | |
| // array index: ADJ | (UP << 2) | (NUNIT << 3) | |
| (Eigen::internal::band_solve_triangular_selector<EIGEN_BLAS_INT, Lower | 0, Scalar, Conj, Scalar, RowMajor>::run), | |
| 0, | |
| // array index: NOTR | (LO << 2) | (NUNIT << 3) | |
| (Eigen::internal::band_solve_triangular_selector<EIGEN_BLAS_INT, Lower | 0, Scalar, false, Scalar, | |
| ColMajor>::run), | |
| // array index: TR | (LO << 2) | (NUNIT << 3) | |
| (Eigen::internal::band_solve_triangular_selector<EIGEN_BLAS_INT, Upper | 0, Scalar, false, Scalar, | |
| RowMajor>::run), | |
| // array index: ADJ | (LO << 2) | (NUNIT << 3) | |
| (Eigen::internal::band_solve_triangular_selector<EIGEN_BLAS_INT, Upper | 0, Scalar, Conj, Scalar, RowMajor>::run), | |
| 0, | |
| // array index: NOTR | (UP << 2) | (UNIT << 3) | |
| (Eigen::internal::band_solve_triangular_selector<EIGEN_BLAS_INT, Upper | UnitDiag, Scalar, false, Scalar, | |
| ColMajor>::run), | |
| // array index: TR | (UP << 2) | (UNIT << 3) | |
| (Eigen::internal::band_solve_triangular_selector<EIGEN_BLAS_INT, Lower | UnitDiag, Scalar, false, Scalar, | |
| RowMajor>::run), | |
| // array index: ADJ | (UP << 2) | (UNIT << 3) | |
| (Eigen::internal::band_solve_triangular_selector<EIGEN_BLAS_INT, Lower | UnitDiag, Scalar, Conj, Scalar, | |
| RowMajor>::run), | |
| 0, | |
| // array index: NOTR | (LO << 2) | (UNIT << 3) | |
| (Eigen::internal::band_solve_triangular_selector<EIGEN_BLAS_INT, Lower | UnitDiag, Scalar, false, Scalar, | |
| ColMajor>::run), | |
| // array index: TR | (LO << 2) | (UNIT << 3) | |
| (Eigen::internal::band_solve_triangular_selector<EIGEN_BLAS_INT, Upper | UnitDiag, Scalar, false, Scalar, | |
| RowMajor>::run), | |
| // array index: ADJ | (LO << 2) | (UNIT << 3) | |
| (Eigen::internal::band_solve_triangular_selector<EIGEN_BLAS_INT, Upper | UnitDiag, Scalar, Conj, Scalar, | |
| RowMajor>::run), | |
| 0, | |
| }; | |
| Scalar *a = reinterpret_cast<Scalar *>(pa); | |
| Scalar *x = reinterpret_cast<Scalar *>(px); | |
| EIGEN_BLAS_INT coeff_rows = *k + 1; | |
| EIGEN_BLAS_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, kBlasNameLength); | |
| if (*n == 0 || (*k == 0 && DIAG(*diag) == UNIT)) return; | |
| EIGEN_BLAS_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, EIGEN_BLAS_INT *n, RealScalar *pap, RealScalar *px, EIGEN_BLAS_INT *incx) { | |
| typedef void (*functype)(EIGEN_BLAS_INT, const Scalar *, const Scalar *, Scalar *, Scalar); | |
| using Eigen::ColMajor; | |
| using Eigen::Lower; | |
| using Eigen::OnTheLeft; | |
| using Eigen::RowMajor; | |
| using Eigen::UnitDiag; | |
| using Eigen::Upper; | |
| static const functype func[16] = { | |
| // array index: NOTR | (UP << 2) | (NUNIT << 3) | |
| (Eigen::internal::packed_triangular_matrix_vector_product<EIGEN_BLAS_INT, Upper | 0, Scalar, false, Scalar, false, | |
| ColMajor>::run), | |
| // array index: TR | (UP << 2) | (NUNIT << 3) | |
| (Eigen::internal::packed_triangular_matrix_vector_product<EIGEN_BLAS_INT, Lower | 0, Scalar, false, Scalar, false, | |
| RowMajor>::run), | |
| // array index: ADJ | (UP << 2) | (NUNIT << 3) | |
| (Eigen::internal::packed_triangular_matrix_vector_product<EIGEN_BLAS_INT, Lower | 0, Scalar, Conj, Scalar, false, | |
| RowMajor>::run), | |
| 0, | |
| // array index: NOTR | (LO << 2) | (NUNIT << 3) | |
| (Eigen::internal::packed_triangular_matrix_vector_product<EIGEN_BLAS_INT, Lower | 0, Scalar, false, Scalar, false, | |
| ColMajor>::run), | |
| // array index: TR | (LO << 2) | (NUNIT << 3) | |
| (Eigen::internal::packed_triangular_matrix_vector_product<EIGEN_BLAS_INT, Upper | 0, Scalar, false, Scalar, false, | |
| RowMajor>::run), | |
| // array index: ADJ | (LO << 2) | (NUNIT << 3) | |
| (Eigen::internal::packed_triangular_matrix_vector_product<EIGEN_BLAS_INT, Upper | 0, Scalar, Conj, Scalar, false, | |
| RowMajor>::run), | |
| 0, | |
| // array index: NOTR | (UP << 2) | (UNIT << 3) | |
| (Eigen::internal::packed_triangular_matrix_vector_product<EIGEN_BLAS_INT, Upper | UnitDiag, Scalar, false, Scalar, | |
| false, ColMajor>::run), | |
| // array index: TR | (UP << 2) | (UNIT << 3) | |
| (Eigen::internal::packed_triangular_matrix_vector_product<EIGEN_BLAS_INT, Lower | UnitDiag, Scalar, false, Scalar, | |
| false, RowMajor>::run), | |
| // array index: ADJ | (UP << 2) | (UNIT << 3) | |
| (Eigen::internal::packed_triangular_matrix_vector_product<EIGEN_BLAS_INT, Lower | UnitDiag, Scalar, Conj, Scalar, | |
| false, RowMajor>::run), | |
| 0, | |
| // array index: NOTR | (LO << 2) | (UNIT << 3) | |
| (Eigen::internal::packed_triangular_matrix_vector_product<EIGEN_BLAS_INT, Lower | UnitDiag, Scalar, false, Scalar, | |
| false, ColMajor>::run), | |
| // array index: TR | (LO << 2) | (UNIT << 3) | |
| (Eigen::internal::packed_triangular_matrix_vector_product<EIGEN_BLAS_INT, Upper | UnitDiag, Scalar, false, Scalar, | |
| false, RowMajor>::run), | |
| // array index: ADJ | (LO << 2) | (UNIT << 3) | |
| (Eigen::internal::packed_triangular_matrix_vector_product<EIGEN_BLAS_INT, Upper | UnitDiag, Scalar, Conj, Scalar, | |
| false, RowMajor>::run), | |
| 0}; | |
| Scalar *ap = reinterpret_cast<Scalar *>(pap); | |
| Scalar *x = reinterpret_cast<Scalar *>(px); | |
| EIGEN_BLAS_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, kBlasNameLength); | |
| if (*n == 0) return; | |
| Scalar *actual_x = get_compact_vector(x, *n, *incx); | |
| Eigen::Matrix<Scalar, Eigen::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, EIGEN_BLAS_INT *n, RealScalar *pap, RealScalar *px, EIGEN_BLAS_INT *incx) { | |
| typedef void (*functype)(EIGEN_BLAS_INT, const Scalar *, Scalar *); | |
| using Eigen::ColMajor; | |
| using Eigen::Lower; | |
| using Eigen::OnTheLeft; | |
| using Eigen::RowMajor; | |
| using Eigen::UnitDiag; | |
| using Eigen::Upper; | |
| static const functype func[16] = { | |
| // array index: NOTR | (UP << 2) | (NUNIT << 3) | |
| (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, EIGEN_BLAS_INT, OnTheLeft, Upper | 0, false, | |
| ColMajor>::run), | |
| // array index: TR | (UP << 2) | (NUNIT << 3) | |
| (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, EIGEN_BLAS_INT, OnTheLeft, Lower | 0, false, | |
| RowMajor>::run), | |
| // array index: ADJ | (UP << 2) | (NUNIT << 3) | |
| (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, EIGEN_BLAS_INT, OnTheLeft, Lower | 0, Conj, | |
| RowMajor>::run), | |
| 0, | |
| // array index: NOTR | (LO << 2) | (NUNIT << 3) | |
| (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, EIGEN_BLAS_INT, OnTheLeft, Lower | 0, false, | |
| ColMajor>::run), | |
| // array index: TR | (LO << 2) | (NUNIT << 3) | |
| (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, EIGEN_BLAS_INT, OnTheLeft, Upper | 0, false, | |
| RowMajor>::run), | |
| // array index: ADJ | (LO << 2) | (NUNIT << 3) | |
| (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, EIGEN_BLAS_INT, OnTheLeft, Upper | 0, Conj, | |
| RowMajor>::run), | |
| 0, | |
| // array index: NOTR | (UP << 2) | (UNIT << 3) | |
| (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, EIGEN_BLAS_INT, OnTheLeft, Upper | UnitDiag, | |
| false, ColMajor>::run), | |
| // array index: TR | (UP << 2) | (UNIT << 3) | |
| (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, EIGEN_BLAS_INT, OnTheLeft, Lower | UnitDiag, | |
| false, RowMajor>::run), | |
| // array index: ADJ | (UP << 2) | (UNIT << 3) | |
| (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, EIGEN_BLAS_INT, OnTheLeft, Lower | UnitDiag, | |
| Conj, RowMajor>::run), | |
| 0, | |
| // array index: NOTR | (LO << 2) | (UNIT << 3) | |
| (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, EIGEN_BLAS_INT, OnTheLeft, Lower | UnitDiag, | |
| false, ColMajor>::run), | |
| // array index: TR | (LO << 2) | (UNIT << 3) | |
| (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, EIGEN_BLAS_INT, OnTheLeft, Upper | UnitDiag, | |
| false, RowMajor>::run), | |
| // array index: ADJ | (LO << 2) | (UNIT << 3) | |
| (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, EIGEN_BLAS_INT, OnTheLeft, Upper | UnitDiag, | |
| Conj, RowMajor>::run), | |
| 0}; | |
| Scalar *ap = reinterpret_cast<Scalar *>(pap); | |
| Scalar *x = reinterpret_cast<Scalar *>(px); | |
| EIGEN_BLAS_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, kBlasNameLength); | |
| 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); | |
| } |