CFD/eigen
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
master
eigen/blas/level2_impl.h

571 lines
25 KiB
C++
Raw Permalink Blame History

// 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;
}
Reference in New Issue View Git Blame Copy Permalink