From 8a73c6490f0a144be2078795537007988e785676 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Antonio=20S=C3=A1nchez?= Date: Thu, 22 Feb 2024 22:51:42 +0000 Subject: [PATCH] Remove "using namespace Eigen" from blas/common.h. --- blas/BandTriangularSolver.h | 4 +- blas/GeneralRank1Update.h | 4 +- blas/PackedSelfadjointProduct.h | 4 +- blas/PackedTriangularMatrixVector.h | 4 +- blas/PackedTriangularSolverVector.h | 4 +- blas/Rank2Update.h | 4 +- blas/common.h | 67 ++--- blas/double.cpp | 2 +- blas/level1_cplx_impl.h | 16 +- blas/level1_impl.h | 4 +- blas/level1_real_impl.h | 14 +- blas/level2_cplx_impl.h | 28 +- blas/level2_impl.h | 226 ++++++++++------ blas/level2_real_impl.h | 29 ++- blas/level3_impl.h | 391 +++++++++++++++++----------- lapack/cholesky.inc | 12 +- lapack/eigenvalues.inc | 5 +- lapack/lu.inc | 2 + lapack/svd.inc | 18 +- 19 files changed, 503 insertions(+), 335 deletions(-) diff --git a/blas/BandTriangularSolver.h b/blas/BandTriangularSolver.h index 014af24e4..f9bfdc1d4 100644 --- a/blas/BandTriangularSolver.h +++ b/blas/BandTriangularSolver.h @@ -10,6 +10,7 @@ #ifndef EIGEN_BAND_TRIANGULARSOLVER_H #define EIGEN_BAND_TRIANGULARSOLVER_H +namespace Eigen { namespace internal { /* \internal @@ -77,6 +78,7 @@ struct band_solve_triangular_selector { } }; -} // end namespace internal +} // namespace internal +} // namespace Eigen #endif // EIGEN_GENERAL_RANK1UPDATE_H diff --git a/blas/PackedSelfadjointProduct.h b/blas/PackedSelfadjointProduct.h index 655da513f..5109960fb 100644 --- a/blas/PackedSelfadjointProduct.h +++ b/blas/PackedSelfadjointProduct.h @@ -10,6 +10,7 @@ #ifndef EIGEN_SELFADJOINT_PACKED_PRODUCT_H #define EIGEN_SELFADJOINT_PACKED_PRODUCT_H +namespace Eigen { namespace internal { /* Optimized matrix += alpha * uv' @@ -45,6 +46,7 @@ struct selfadjoint_packed_rank1_update @@ -69,6 +70,7 @@ struct packed_triangular_solve_vector::Real RealScalar; +typedef Eigen::NumTraits::Real RealScalar; typedef std::complex Complex; enum { IsComplex = Eigen::NumTraits::IsComplex, Conj = IsComplex }; -typedef Matrix PlainMatrixType; -typedef Map, 0, OuterStride<> > MatrixType; -typedef Map, 0, OuterStride<> > ConstMatrixType; -typedef Map, 0, InnerStride > StridedVectorType; -typedef Map > CompactVectorType; +typedef Eigen::Matrix PlainMatrixType; +typedef Eigen::Map, 0, Eigen::OuterStride<> > + MatrixType; +typedef Eigen::Map, 0, + Eigen::OuterStride<> > + ConstMatrixType; +typedef Eigen::Map, 0, Eigen::InnerStride > StridedVectorType; +typedef Eigen::Map > CompactVectorType; template -Map, 0, OuterStride<> > matrix(T* data, int rows, int cols, int stride) { - return Map, 0, OuterStride<> >(data, rows, cols, OuterStride<>(stride)); +Eigen::Map, 0, Eigen::OuterStride<> > matrix( + T* data, int rows, int cols, int stride) { + return Eigen::Map, 0, Eigen::OuterStride<> >( + data, rows, cols, Eigen::OuterStride<>(stride)); } template -Map, 0, OuterStride<> > matrix(const T* data, int rows, int cols, - int stride) { - return Map, 0, OuterStride<> >(data, rows, cols, OuterStride<>(stride)); +Eigen::Map, 0, Eigen::OuterStride<> > matrix( + const T* data, int rows, int cols, int stride) { + return Eigen::Map, 0, Eigen::OuterStride<> >( + data, rows, cols, Eigen::OuterStride<>(stride)); } template -Map, 0, InnerStride > make_vector(T* data, int size, int incr) { - return Map, 0, InnerStride >(data, size, InnerStride(incr)); +Eigen::Map, 0, Eigen::InnerStride > make_vector(T* data, int size, + int incr) { + return Eigen::Map, 0, Eigen::InnerStride >( + data, size, Eigen::InnerStride(incr)); } template -Map, 0, InnerStride > make_vector(const T* data, int size, int incr) { - return Map, 0, InnerStride >(data, size, InnerStride(incr)); +Eigen::Map, 0, Eigen::InnerStride > make_vector(const T* data, + int size, + int incr) { + return Eigen::Map, 0, Eigen::InnerStride >( + data, size, Eigen::InnerStride(incr)); } template -Map > make_vector(T* data, int size) { - return Map >(data, size); +Eigen::Map > make_vector(T* data, int size) { + return Eigen::Map >(data, size); } template -Map > make_vector(const T* data, int size) { - return Map >(data, size); +Eigen::Map > make_vector(const T* data, int size) { + return Eigen::Map >(data, size); } template diff --git a/blas/double.cpp b/blas/double.cpp index 429866566..28a25634b 100644 --- a/blas/double.cpp +++ b/blas/double.cpp @@ -19,7 +19,7 @@ #include "level2_real_impl.h" #include "level3_impl.h" -double EIGEN_BLAS_FUNC_NAME(sdot)(int* n, float* x, int* incx, float* y, int* incy) { +extern "C" double EIGEN_BLAS_FUNC_NAME(sdot)(int* n, float* x, int* incx, float* y, int* incy) { if (*n <= 0) return 0; if (*incx == 1 && *incy == 1) diff --git a/blas/level1_cplx_impl.h b/blas/level1_cplx_impl.h index be88b92b1..3181a5038 100644 --- a/blas/level1_cplx_impl.h +++ b/blas/level1_cplx_impl.h @@ -11,7 +11,7 @@ struct scalar_norm1_op { typedef RealScalar result_type; - inline RealScalar operator()(const Scalar &a) const { return numext::norm1(a); } + inline RealScalar operator()(const Scalar &a) const { return Eigen::numext::norm1(a); } }; namespace Eigen { namespace internal { @@ -40,7 +40,7 @@ extern "C" int EIGEN_CAT(i, EIGEN_BLAS_FUNC_NAME(amax))(int *n, RealScalar *px, if (*n <= 0) return 0; Scalar *x = reinterpret_cast(px); - DenseIndex ret; + Eigen::DenseIndex ret; if (*incx == 1) make_vector(x, *n).unaryExpr().maxCoeff(&ret); else @@ -52,7 +52,7 @@ extern "C" int EIGEN_CAT(i, EIGEN_BLAS_FUNC_NAME(amin))(int *n, RealScalar *px, if (*n <= 0) return 0; Scalar *x = reinterpret_cast(px); - DenseIndex ret; + Eigen::DenseIndex ret; if (*incx == 1) make_vector(x, *n).unaryExpr().minCoeff(&ret); else @@ -132,16 +132,16 @@ EIGEN_BLAS_FUNC(EIGEN_CAT(REAL_SCALAR_SUFFIX, rot)) StridedVectorType vx(make_vector(x, *n, std::abs(*incx))); StridedVectorType vy(make_vector(y, *n, std::abs(*incy))); - Reverse rvx(vx); - Reverse rvy(vy); + Eigen::Reverse rvx(vx); + Eigen::Reverse rvy(vy); // TODO implement mixed real-scalar rotations if (*incx < 0 && *incy > 0) - internal::apply_rotation_in_the_plane(rvx, vy, JacobiRotation(c, s)); + Eigen::internal::apply_rotation_in_the_plane(rvx, vy, Eigen::JacobiRotation(c, s)); else if (*incx > 0 && *incy < 0) - internal::apply_rotation_in_the_plane(vx, rvy, JacobiRotation(c, s)); + Eigen::internal::apply_rotation_in_the_plane(vx, rvy, Eigen::JacobiRotation(c, s)); else - internal::apply_rotation_in_the_plane(vx, vy, JacobiRotation(c, s)); + Eigen::internal::apply_rotation_in_the_plane(vx, vy, Eigen::JacobiRotation(c, s)); } EIGEN_BLAS_FUNC(EIGEN_CAT(REAL_SCALAR_SUFFIX, scal))(int *n, RealScalar *palpha, RealScalar *px, int *incx) { diff --git a/blas/level1_impl.h b/blas/level1_impl.h index 2422d1017..a65af92d5 100644 --- a/blas/level1_impl.h +++ b/blas/level1_impl.h @@ -88,10 +88,10 @@ EIGEN_BLAS_FUNC(rotg)(RealScalar *pa, RealScalar *pb, RealScalar *pc, RealScalar a = b; } else { scale = abs(a) + abs(b); - norm = scale * sqrt((numext::abs2(a / scale)) + (numext::abs2(b / scale))); + norm = scale * sqrt((Eigen::numext::abs2(a / scale)) + (Eigen::numext::abs2(b / scale))); alpha = a / abs(a); *c = abs(a) / norm; - *s = alpha * numext::conj(b) / norm; + *s = alpha * Eigen::numext::conj(b) / norm; a = alpha * norm; } #endif diff --git a/blas/level1_real_impl.h b/blas/level1_real_impl.h index cd9c1899b..202f4324d 100644 --- a/blas/level1_real_impl.h +++ b/blas/level1_real_impl.h @@ -28,7 +28,7 @@ extern "C" int EIGEN_CAT(i, EIGEN_BLAS_FUNC_NAME(amax))(int *n, Scalar *px, int if (*n <= 0) return 0; Scalar *x = reinterpret_cast(px); - DenseIndex ret; + Eigen::DenseIndex ret; if (*incx == 1) make_vector(x, *n).cwiseAbs().maxCoeff(&ret); else @@ -40,7 +40,7 @@ extern "C" int EIGEN_CAT(i, EIGEN_BLAS_FUNC_NAME(amin))(int *n, Scalar *px, int if (*n <= 0) return 0; Scalar *x = reinterpret_cast(px); - DenseIndex ret; + Eigen::DenseIndex ret; if (*incx == 1) make_vector(x, *n).cwiseAbs().minCoeff(&ret); else @@ -97,15 +97,15 @@ EIGEN_BLAS_FUNC(rot)(int *n, Scalar *px, int *incx, Scalar *py, int *incy, Scala StridedVectorType vx(make_vector(x, *n, std::abs(*incx))); StridedVectorType vy(make_vector(y, *n, std::abs(*incy))); - Reverse rvx(vx); - Reverse rvy(vy); + Eigen::Reverse rvx(vx); + Eigen::Reverse rvy(vy); if (*incx < 0 && *incy > 0) - internal::apply_rotation_in_the_plane(rvx, vy, JacobiRotation(c, s)); + Eigen::internal::apply_rotation_in_the_plane(rvx, vy, Eigen::JacobiRotation(c, s)); else if (*incx > 0 && *incy < 0) - internal::apply_rotation_in_the_plane(vx, rvy, JacobiRotation(c, s)); + Eigen::internal::apply_rotation_in_the_plane(vx, rvy, Eigen::JacobiRotation(c, s)); else - internal::apply_rotation_in_the_plane(vx, vy, JacobiRotation(c, s)); + Eigen::internal::apply_rotation_in_the_plane(vx, vy, Eigen::JacobiRotation(c, s)); } /* diff --git a/blas/level2_cplx_impl.h b/blas/level2_cplx_impl.h index f04dda120..d1ce492ac 100644 --- a/blas/level2_cplx_impl.h +++ b/blas/level2_cplx_impl.h @@ -22,9 +22,11 @@ EIGEN_BLAS_FUNC(hemv) typedef void (*functype)(int, const Scalar *, int, const Scalar *, Scalar *, Scalar); static const functype func[2] = { // array index: UP - (internal::selfadjoint_matrix_vector_product::run), + (Eigen::internal::selfadjoint_matrix_vector_product::run), // array index: LO - (internal::selfadjoint_matrix_vector_product::run), + (Eigen::internal::selfadjoint_matrix_vector_product::run), }; const Scalar *a = reinterpret_cast(pa); @@ -107,9 +109,9 @@ EIGEN_BLAS_FUNC(hpr)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int typedef void (*functype)(int, Scalar *, const Scalar *, RealScalar); static const functype func[2] = { // array index: UP - (internal::selfadjoint_packed_rank1_update::run), + (Eigen::internal::selfadjoint_packed_rank1_update::run), // array index: LO - (internal::selfadjoint_packed_rank1_update::run), + (Eigen::internal::selfadjoint_packed_rank1_update::run), }; Scalar *x = reinterpret_cast(px); @@ -149,9 +151,9 @@ EIGEN_BLAS_FUNC(hpr2) typedef void (*functype)(int, Scalar *, const Scalar *, const Scalar *, Scalar); static const functype func[2] = { // array index: UP - (internal::packed_rank2_update_selector::run), + (Eigen::internal::packed_rank2_update_selector::run), // array index: LO - (internal::packed_rank2_update_selector::run), + (Eigen::internal::packed_rank2_update_selector::run), }; Scalar *x = reinterpret_cast(px); @@ -195,9 +197,9 @@ EIGEN_BLAS_FUNC(her)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int typedef void (*functype)(int, Scalar *, int, const Scalar *, const Scalar *, const Scalar &); static const functype func[2] = { // array index: UP - (selfadjoint_rank1_update::run), + (Eigen::selfadjoint_rank1_update::run), // array index: LO - (selfadjoint_rank1_update::run), + (Eigen::selfadjoint_rank1_update::run), }; Scalar *x = reinterpret_cast(px); @@ -242,9 +244,9 @@ EIGEN_BLAS_FUNC(her2) typedef void (*functype)(int, Scalar *, int, const Scalar *, const Scalar *, Scalar); static const functype func[2] = { // array index: UP - (internal::rank2_update_selector::run), + (Eigen::internal::rank2_update_selector::run), // array index: LO - (internal::rank2_update_selector::run), + (Eigen::internal::rank2_update_selector::run), }; Scalar *x = reinterpret_cast(px); @@ -313,7 +315,8 @@ EIGEN_BLAS_FUNC(geru) Scalar *x_cpy = get_compact_vector(x, *m, *incx); Scalar *y_cpy = get_compact_vector(y, *n, *incy); - internal::general_rank1_update::run(*m, *n, a, *lda, x_cpy, y_cpy, alpha); + Eigen::internal::general_rank1_update::run(*m, *n, a, *lda, x_cpy, y_cpy, + alpha); if (x_cpy != x) delete[] x_cpy; if (y_cpy != y) delete[] y_cpy; @@ -351,7 +354,8 @@ EIGEN_BLAS_FUNC(gerc) Scalar *x_cpy = get_compact_vector(x, *m, *incx); Scalar *y_cpy = get_compact_vector(y, *n, *incy); - internal::general_rank1_update::run(*m, *n, a, *lda, x_cpy, y_cpy, alpha); + Eigen::internal::general_rank1_update::run(*m, *n, a, *lda, x_cpy, y_cpy, + alpha); if (x_cpy != x) delete[] x_cpy; if (y_cpy != y) delete[] y_cpy; diff --git a/blas/level2_impl.h b/blas/level2_impl.h index 5721ee603..ca9f48f3d 100644 --- a/blas/level2_impl.h +++ b/blas/level2_impl.h @@ -13,12 +13,13 @@ template LhsMapper; - typedef internal::const_blas_data_mapper RhsMapper; + typedef Eigen::internal::const_blas_data_mapper LhsMapper; + typedef Eigen::internal::const_blas_data_mapper RhsMapper; - internal::general_matrix_vector_product::run(rows, cols, LhsMapper(lhs, lhsStride), - RhsMapper(rhs, rhsIncr), res, resIncr, alpha); + Eigen::internal::general_matrix_vector_product::run(rows, cols, LhsMapper(lhs, lhsStride), + RhsMapper(rhs, rhsIncr), res, resIncr, + alpha); } }; @@ -26,12 +27,13 @@ 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::run), - // array index: TR - (general_matrix_vector_product_wrapper::run), - // array index: ADJ - (general_matrix_vector_product_wrapper::run), 0}; + static const functype func[4] = { + // array index: NOTR + (general_matrix_vector_product_wrapper::run), + // array index: TR + (general_matrix_vector_product_wrapper::run), + // array index: ADJ + (general_matrix_vector_product_wrapper::run), 0}; const Scalar *a = reinterpret_cast(pa); const Scalar *b = reinterpret_cast(pb); @@ -84,31 +86,43 @@ 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 *); + 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) - (internal::triangular_solve_vector::run), + (Eigen::internal::triangular_solve_vector::run), // array index: TR | (UP << 2) | (NUNIT << 3) - (internal::triangular_solve_vector::run), + (Eigen::internal::triangular_solve_vector::run), // array index: ADJ | (UP << 2) | (NUNIT << 3) - (internal::triangular_solve_vector::run), 0, + (Eigen::internal::triangular_solve_vector::run), 0, // array index: NOTR | (LO << 2) | (NUNIT << 3) - (internal::triangular_solve_vector::run), + (Eigen::internal::triangular_solve_vector::run), // array index: TR | (LO << 2) | (NUNIT << 3) - (internal::triangular_solve_vector::run), + (Eigen::internal::triangular_solve_vector::run), // array index: ADJ | (LO << 2) | (NUNIT << 3) - (internal::triangular_solve_vector::run), 0, + (Eigen::internal::triangular_solve_vector::run), 0, // array index: NOTR | (UP << 2) | (UNIT << 3) - (internal::triangular_solve_vector::run), + (Eigen::internal::triangular_solve_vector::run), // array index: TR | (UP << 2) | (UNIT << 3) - (internal::triangular_solve_vector::run), + (Eigen::internal::triangular_solve_vector::run), // array index: ADJ | (UP << 2) | (UNIT << 3) - (internal::triangular_solve_vector::run), 0, + (Eigen::internal::triangular_solve_vector::run), + 0, // array index: NOTR | (LO << 2) | (UNIT << 3) - (internal::triangular_solve_vector::run), + (Eigen::internal::triangular_solve_vector::run), // array index: TR | (LO << 2) | (UNIT << 3) - (internal::triangular_solve_vector::run), + (Eigen::internal::triangular_solve_vector::run), // array index: ADJ | (LO << 2) | (UNIT << 3) - (internal::triangular_solve_vector::run), 0}; + (Eigen::internal::triangular_solve_vector::run), + 0}; const Scalar *a = reinterpret_cast(pa); Scalar *b = reinterpret_cast(pb); @@ -140,32 +154,46 @@ 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 &); + 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) - (internal::triangular_matrix_vector_product::run), + (Eigen::internal::triangular_matrix_vector_product::run), // array index: TR | (UP << 2) | (NUNIT << 3) - (internal::triangular_matrix_vector_product::run), + (Eigen::internal::triangular_matrix_vector_product::run), // array index: ADJ | (UP << 2) | (NUNIT << 3) - (internal::triangular_matrix_vector_product::run), 0, + (Eigen::internal::triangular_matrix_vector_product::run), + 0, // array index: NOTR | (LO << 2) | (NUNIT << 3) - (internal::triangular_matrix_vector_product::run), + (Eigen::internal::triangular_matrix_vector_product::run), // array index: TR | (LO << 2) | (NUNIT << 3) - (internal::triangular_matrix_vector_product::run), + (Eigen::internal::triangular_matrix_vector_product::run), // array index: ADJ | (LO << 2) | (NUNIT << 3) - (internal::triangular_matrix_vector_product::run), 0, + (Eigen::internal::triangular_matrix_vector_product::run), + 0, // array index: NOTR | (UP << 2) | (UNIT << 3) - (internal::triangular_matrix_vector_product::run), + (Eigen::internal::triangular_matrix_vector_product::run), // array index: TR | (UP << 2) | (UNIT << 3) - (internal::triangular_matrix_vector_product::run), + (Eigen::internal::triangular_matrix_vector_product::run), // array index: ADJ | (UP << 2) | (UNIT << 3) - (internal::triangular_matrix_vector_product::run), + (Eigen::internal::triangular_matrix_vector_product::run), 0, // array index: NOTR | (LO << 2) | (UNIT << 3) - (internal::triangular_matrix_vector_product::run), + (Eigen::internal::triangular_matrix_vector_product::run), // array index: TR | (LO << 2) | (UNIT << 3) - (internal::triangular_matrix_vector_product::run), + (Eigen::internal::triangular_matrix_vector_product::run), // array index: ADJ | (LO << 2) | (UNIT << 3) - (internal::triangular_matrix_vector_product::run), + (Eigen::internal::triangular_matrix_vector_product::run), 0}; const Scalar *a = reinterpret_cast(pa); @@ -189,7 +217,7 @@ EIGEN_BLAS_FUNC(trmv) if (*n == 0) return; Scalar *actual_b = get_compact_vector(b, *n, *incb); - Matrix res(*n); + Eigen::Matrix res(*n); res.setZero(); int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3); @@ -345,34 +373,40 @@ EIGEN_BLAS_FUNC(tbmv)(char *uplo, char *opa, char *diag, int *n, int *k, RealSca 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 *); + 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) - (internal::band_solve_triangular_selector::run), + (Eigen::internal::band_solve_triangular_selector::run), // array index: TR | (UP << 2) | (NUNIT << 3) - (internal::band_solve_triangular_selector::run), + (Eigen::internal::band_solve_triangular_selector::run), // array index: ADJ | (UP << 2) | (NUNIT << 3) - (internal::band_solve_triangular_selector::run), + (Eigen::internal::band_solve_triangular_selector::run), 0, // array index: NOTR | (LO << 2) | (NUNIT << 3) - (internal::band_solve_triangular_selector::run), + (Eigen::internal::band_solve_triangular_selector::run), // array index: TR | (LO << 2) | (NUNIT << 3) - (internal::band_solve_triangular_selector::run), + (Eigen::internal::band_solve_triangular_selector::run), // array index: ADJ | (LO << 2) | (NUNIT << 3) - (internal::band_solve_triangular_selector::run), + (Eigen::internal::band_solve_triangular_selector::run), 0, // array index: NOTR | (UP << 2) | (UNIT << 3) - (internal::band_solve_triangular_selector::run), + (Eigen::internal::band_solve_triangular_selector::run), // array index: TR | (UP << 2) | (UNIT << 3) - (internal::band_solve_triangular_selector::run), + (Eigen::internal::band_solve_triangular_selector::run), // array index: ADJ | (UP << 2) | (UNIT << 3) - (internal::band_solve_triangular_selector::run), + (Eigen::internal::band_solve_triangular_selector::run), 0, // array index: NOTR | (LO << 2) | (UNIT << 3) - (internal::band_solve_triangular_selector::run), + (Eigen::internal::band_solve_triangular_selector::run), // array index: TR | (LO << 2) | (UNIT << 3) - (internal::band_solve_triangular_selector::run), + (Eigen::internal::band_solve_triangular_selector::run), // array index: ADJ | (LO << 2) | (UNIT << 3) - (internal::band_solve_triangular_selector::run), + (Eigen::internal::band_solve_triangular_selector::run), 0, }; @@ -420,40 +454,52 @@ EIGEN_BLAS_FUNC(tbsv) */ 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); + 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) - (internal::packed_triangular_matrix_vector_product::run), + (Eigen::internal::packed_triangular_matrix_vector_product::run), // array index: TR | (UP << 2) | (NUNIT << 3) - (internal::packed_triangular_matrix_vector_product::run), + (Eigen::internal::packed_triangular_matrix_vector_product::run), // array index: ADJ | (UP << 2) | (NUNIT << 3) - (internal::packed_triangular_matrix_vector_product::run), + (Eigen::internal::packed_triangular_matrix_vector_product::run), 0, // array index: NOTR | (LO << 2) | (NUNIT << 3) - (internal::packed_triangular_matrix_vector_product::run), + (Eigen::internal::packed_triangular_matrix_vector_product::run), // array index: TR | (LO << 2) | (NUNIT << 3) - (internal::packed_triangular_matrix_vector_product::run), + (Eigen::internal::packed_triangular_matrix_vector_product::run), // array index: ADJ | (LO << 2) | (NUNIT << 3) - (internal::packed_triangular_matrix_vector_product::run), + (Eigen::internal::packed_triangular_matrix_vector_product::run), 0, // array index: NOTR | (UP << 2) | (UNIT << 3) - (internal::packed_triangular_matrix_vector_product::run), + (Eigen::internal::packed_triangular_matrix_vector_product::run), // array index: TR | (UP << 2) | (UNIT << 3) - (internal::packed_triangular_matrix_vector_product::run), + (Eigen::internal::packed_triangular_matrix_vector_product::run), // array index: ADJ | (UP << 2) | (UNIT << 3) - (internal::packed_triangular_matrix_vector_product::run), + (Eigen::internal::packed_triangular_matrix_vector_product::run), 0, // array index: NOTR | (LO << 2) | (UNIT << 3) - (internal::packed_triangular_matrix_vector_product::run), + (Eigen::internal::packed_triangular_matrix_vector_product::run), // array index: TR | (LO << 2) | (UNIT << 3) - (internal::packed_triangular_matrix_vector_product::run), + (Eigen::internal::packed_triangular_matrix_vector_product::run), // array index: ADJ | (LO << 2) | (UNIT << 3) - (internal::packed_triangular_matrix_vector_product::run), + (Eigen::internal::packed_triangular_matrix_vector_product::run), 0}; Scalar *ap = reinterpret_cast(pap); @@ -475,7 +521,7 @@ EIGEN_BLAS_FUNC(tpmv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pap if (*n == 0) return; Scalar *actual_x = get_compact_vector(x, *n, *incx); - Matrix res(*n); + Eigen::Matrix res(*n); res.setZero(); int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3); @@ -499,36 +545,50 @@ EIGEN_BLAS_FUNC(tpmv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pap */ 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 *); + 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) - (internal::packed_triangular_solve_vector::run), + (Eigen::internal::packed_triangular_solve_vector::run), // array index: TR | (UP << 2) | (NUNIT << 3) - (internal::packed_triangular_solve_vector::run), + (Eigen::internal::packed_triangular_solve_vector::run), // array index: ADJ | (UP << 2) | (NUNIT << 3) - (internal::packed_triangular_solve_vector::run), 0, + (Eigen::internal::packed_triangular_solve_vector::run), + 0, // array index: NOTR | (LO << 2) | (NUNIT << 3) - (internal::packed_triangular_solve_vector::run), + (Eigen::internal::packed_triangular_solve_vector::run), // array index: TR | (LO << 2) | (NUNIT << 3) - (internal::packed_triangular_solve_vector::run), + (Eigen::internal::packed_triangular_solve_vector::run), // array index: ADJ | (LO << 2) | (NUNIT << 3) - (internal::packed_triangular_solve_vector::run), 0, + (Eigen::internal::packed_triangular_solve_vector::run), + 0, // array index: NOTR | (UP << 2) | (UNIT << 3) - (internal::packed_triangular_solve_vector::run), + (Eigen::internal::packed_triangular_solve_vector::run), // array index: TR | (UP << 2) | (UNIT << 3) - (internal::packed_triangular_solve_vector::run), + (Eigen::internal::packed_triangular_solve_vector::run), // array index: ADJ | (UP << 2) | (UNIT << 3) - (internal::packed_triangular_solve_vector::run), + (Eigen::internal::packed_triangular_solve_vector::run), 0, // array index: NOTR | (LO << 2) | (UNIT << 3) - (internal::packed_triangular_solve_vector::run), + (Eigen::internal::packed_triangular_solve_vector::run), // array index: TR | (LO << 2) | (UNIT << 3) - (internal::packed_triangular_solve_vector::run), + (Eigen::internal::packed_triangular_solve_vector::run), // array index: ADJ | (LO << 2) | (UNIT << 3) - (internal::packed_triangular_solve_vector::run), + (Eigen::internal::packed_triangular_solve_vector::run), 0}; Scalar *ap = reinterpret_cast(pap); diff --git a/blas/level2_real_impl.h b/blas/level2_real_impl.h index 5653767d2..415944cc0 100644 --- a/blas/level2_real_impl.h +++ b/blas/level2_real_impl.h @@ -14,11 +14,14 @@ EIGEN_BLAS_FUNC(symv) (const char *uplo, const int *n, const RealScalar *palpha, const RealScalar *pa, const int *lda, const RealScalar *px, const int *incx, const RealScalar *pbeta, RealScalar *py, const int *incy) { typedef void (*functype)(int, const Scalar *, int, const Scalar *, Scalar *, Scalar); + using Eigen::ColMajor; + using Eigen::Lower; + using Eigen::Upper; static const functype func[2] = { // array index: UP - (internal::selfadjoint_matrix_vector_product::run), + (Eigen::internal::selfadjoint_matrix_vector_product::run), // array index: LO - (internal::selfadjoint_matrix_vector_product::run), + (Eigen::internal::selfadjoint_matrix_vector_product::run), }; const Scalar *a = reinterpret_cast(pa); @@ -67,11 +70,14 @@ EIGEN_BLAS_FUNC(syr) (const char *uplo, const int *n, const RealScalar *palpha, const RealScalar *px, const int *incx, RealScalar *pc, const int *ldc) { typedef void (*functype)(int, Scalar *, int, const Scalar *, const Scalar *, const Scalar &); + using Eigen::ColMajor; + using Eigen::Lower; + using Eigen::Upper; static const functype func[2] = { // array index: UP - (selfadjoint_rank1_update::run), + (Eigen::selfadjoint_rank1_update::run), // array index: LO - (selfadjoint_rank1_update::run), + (Eigen::selfadjoint_rank1_update::run), }; const Scalar *x = reinterpret_cast(px); @@ -109,9 +115,9 @@ EIGEN_BLAS_FUNC(syr2) typedef void (*functype)(int, Scalar *, int, const Scalar *, const Scalar *, Scalar); static const functype func[2] = { // array index: UP - (internal::rank2_update_selector::run), + (Eigen::internal::rank2_update_selector::run), // array index: LO - (internal::rank2_update_selector::run), + (Eigen::internal::rank2_update_selector::run), }; const Scalar *x = reinterpret_cast(px); @@ -190,9 +196,9 @@ EIGEN_BLAS_FUNC(spr)(char *uplo, int *n, Scalar *palpha, Scalar *px, int *incx, typedef void (*functype)(int, Scalar *, const Scalar *, Scalar); static const functype func[2] = { // array index: UP - (internal::selfadjoint_packed_rank1_update::run), + (Eigen::internal::selfadjoint_packed_rank1_update::run), // array index: LO - (internal::selfadjoint_packed_rank1_update::run), + (Eigen::internal::selfadjoint_packed_rank1_update::run), }; Scalar *x = reinterpret_cast(px); @@ -232,9 +238,9 @@ EIGEN_BLAS_FUNC(spr2) typedef void (*functype)(int, Scalar *, const Scalar *, const Scalar *, Scalar); static const functype func[2] = { // array index: UP - (internal::packed_rank2_update_selector::run), + (Eigen::internal::packed_rank2_update_selector::run), // array index: LO - (internal::packed_rank2_update_selector::run), + (Eigen::internal::packed_rank2_update_selector::run), }; Scalar *x = reinterpret_cast(px); @@ -299,7 +305,8 @@ EIGEN_BLAS_FUNC(ger) Scalar *x_cpy = get_compact_vector(x, *m, *incx); Scalar *y_cpy = get_compact_vector(y, *n, *incy); - internal::general_rank1_update::run(*m, *n, a, *lda, x_cpy, y_cpy, alpha); + Eigen::internal::general_rank1_update::run(*m, *n, a, *lda, x_cpy, y_cpy, + alpha); if (x_cpy != x) delete[] x_cpy; if (y_cpy != y) delete[] y_cpy; diff --git a/blas/level3_impl.h b/blas/level3_impl.h index a6ddf2614..66a7d468b 100644 --- a/blas/level3_impl.h +++ b/blas/level3_impl.h @@ -15,39 +15,43 @@ EIGEN_BLAS_FUNC(gemm) const int *ldc) { // std::cerr << "in gemm " << *opa << " " << *opb << " " << *m << " " << *n << " " << *k << " " << *lda << " " << // *ldb << " " << *ldc << " " << *palpha << " " << *pbeta << "\n"; + using Eigen::ColMajor; + using Eigen::DenseIndex; + using Eigen::Dynamic; + using Eigen::RowMajor; typedef void (*functype)(DenseIndex, DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, DenseIndex, Scalar, Eigen::internal::level3_blocking &, Eigen::internal::GemmParallelInfo *); static const functype func[12] = { // array index: NOTR | (NOTR << 2) - (internal::general_matrix_matrix_product::run), + (Eigen::internal::general_matrix_matrix_product::run), // array index: TR | (NOTR << 2) - (internal::general_matrix_matrix_product::run), + (Eigen::internal::general_matrix_matrix_product::run), // array index: ADJ | (NOTR << 2) - (internal::general_matrix_matrix_product::run), + (Eigen::internal::general_matrix_matrix_product::run), 0, // array index: NOTR | (TR << 2) - (internal::general_matrix_matrix_product::run), + (Eigen::internal::general_matrix_matrix_product::run), // array index: TR | (TR << 2) - (internal::general_matrix_matrix_product::run), + (Eigen::internal::general_matrix_matrix_product::run), // array index: ADJ | (TR << 2) - (internal::general_matrix_matrix_product::run), + (Eigen::internal::general_matrix_matrix_product::run), 0, // array index: NOTR | (ADJ << 2) - (internal::general_matrix_matrix_product::run), + (Eigen::internal::general_matrix_matrix_product::run), // array index: TR | (ADJ << 2) - (internal::general_matrix_matrix_product::run), + (Eigen::internal::general_matrix_matrix_product::run), // array index: ADJ | (ADJ << 2) - (internal::general_matrix_matrix_product::run), + (Eigen::internal::general_matrix_matrix_product::run), 0}; const Scalar *a = reinterpret_cast(pa); @@ -86,7 +90,8 @@ EIGEN_BLAS_FUNC(gemm) if (*k == 0) return; - internal::gemm_blocking_space blocking(*m, *n, *k, 1, true); + Eigen::internal::gemm_blocking_space blocking(*m, *n, *k, 1, + true); int code = OP(*opa) | (OP(*opb) << 2); func[code](*m, *n, *k, a, *lda, b, *ldb, c, 1, *ldc, alpha, blocking, 0); @@ -97,76 +102,97 @@ EIGEN_BLAS_FUNC(trsm) const RealScalar *palpha, const RealScalar *pa, const int *lda, RealScalar *pb, const int *ldb) { // std::cerr << "in trsm " << *side << " " << *uplo << " " << *opa << " " << *diag << " " << *m << "," << *n << " " // << *palpha << " " << *lda << " " << *ldb<< "\n"; + using Eigen::ColMajor; + using Eigen::DenseIndex; + using Eigen::Dynamic; + using Eigen::Lower; + using Eigen::OnTheLeft; + using Eigen::OnTheRight; + using Eigen::RowMajor; + using Eigen::UnitDiag; + using Eigen::Upper; typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, DenseIndex, Eigen::internal::level3_blocking &); static const functype func[32] = { // array index: NOTR | (LEFT << 2) | (UP << 3) | (NUNIT << 4) - (internal::triangular_solve_matrix::run), + (Eigen::internal::triangular_solve_matrix::run), // array index: TR | (LEFT << 2) | (UP << 3) | (NUNIT << 4) - (internal::triangular_solve_matrix::run), + (Eigen::internal::triangular_solve_matrix::run), // array index: ADJ | (LEFT << 2) | (UP << 3) | (NUNIT << 4) - (internal::triangular_solve_matrix::run), + (Eigen::internal::triangular_solve_matrix::run), 0, // array index: NOTR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4) - (internal::triangular_solve_matrix::run), + (Eigen::internal::triangular_solve_matrix::run), // array index: TR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4) - (internal::triangular_solve_matrix::run), + (Eigen::internal::triangular_solve_matrix::run), // array index: ADJ | (RIGHT << 2) | (UP << 3) | (NUNIT << 4) - (internal::triangular_solve_matrix::run), + (Eigen::internal::triangular_solve_matrix::run), 0, // array index: NOTR | (LEFT << 2) | (LO << 3) | (NUNIT << 4) - (internal::triangular_solve_matrix::run), + (Eigen::internal::triangular_solve_matrix::run), // array index: TR | (LEFT << 2) | (LO << 3) | (NUNIT << 4) - (internal::triangular_solve_matrix::run), + (Eigen::internal::triangular_solve_matrix::run), // array index: ADJ | (LEFT << 2) | (LO << 3) | (NUNIT << 4) - (internal::triangular_solve_matrix::run), + (Eigen::internal::triangular_solve_matrix::run), 0, // array index: NOTR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4) - (internal::triangular_solve_matrix::run), + (Eigen::internal::triangular_solve_matrix::run), // array index: TR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4) - (internal::triangular_solve_matrix::run), + (Eigen::internal::triangular_solve_matrix::run), // array index: ADJ | (RIGHT << 2) | (LO << 3) | (NUNIT << 4) - (internal::triangular_solve_matrix::run), + (Eigen::internal::triangular_solve_matrix::run), 0, // array index: NOTR | (LEFT << 2) | (UP << 3) | (UNIT << 4) - (internal::triangular_solve_matrix::run), + (Eigen::internal::triangular_solve_matrix::run), // array index: TR | (LEFT << 2) | (UP << 3) | (UNIT << 4) - (internal::triangular_solve_matrix::run), + (Eigen::internal::triangular_solve_matrix::run), // array index: ADJ | (LEFT << 2) | (UP << 3) | (UNIT << 4) - (internal::triangular_solve_matrix::run), + (Eigen::internal::triangular_solve_matrix::run), 0, // array index: NOTR | (RIGHT << 2) | (UP << 3) | (UNIT << 4) - (internal::triangular_solve_matrix::run), + (Eigen::internal::triangular_solve_matrix::run), // array index: TR | (RIGHT << 2) | (UP << 3) | (UNIT << 4) - (internal::triangular_solve_matrix::run), + (Eigen::internal::triangular_solve_matrix::run), // array index: ADJ | (RIGHT << 2) | (UP << 3) | (UNIT << 4) - (internal::triangular_solve_matrix::run), + (Eigen::internal::triangular_solve_matrix::run), 0, // array index: NOTR | (LEFT << 2) | (LO << 3) | (UNIT << 4) - (internal::triangular_solve_matrix::run), + (Eigen::internal::triangular_solve_matrix::run), // array index: TR | (LEFT << 2) | (LO << 3) | (UNIT << 4) - (internal::triangular_solve_matrix::run), + (Eigen::internal::triangular_solve_matrix::run), // array index: ADJ | (LEFT << 2) | (LO << 3) | (UNIT << 4) - (internal::triangular_solve_matrix::run), + (Eigen::internal::triangular_solve_matrix::run), 0, // array index: NOTR | (RIGHT << 2) | (LO << 3) | (UNIT << 4) - (internal::triangular_solve_matrix::run), + (Eigen::internal::triangular_solve_matrix::run), // array index: TR | (RIGHT << 2) | (LO << 3) | (UNIT << 4) - (internal::triangular_solve_matrix::run), + (Eigen::internal::triangular_solve_matrix::run), // array index: ADJ | (RIGHT << 2) | (LO << 3) | (UNIT << 4) - (internal::triangular_solve_matrix::run), + (Eigen::internal::triangular_solve_matrix::run), 0}; const Scalar *a = reinterpret_cast(pa); @@ -197,12 +223,12 @@ EIGEN_BLAS_FUNC(trsm) int code = OP(*opa) | (SIDE(*side) << 2) | (UPLO(*uplo) << 3) | (DIAG(*diag) << 4); if (SIDE(*side) == LEFT) { - internal::gemm_blocking_space blocking(*m, *n, *m, 1, - false); + Eigen::internal::gemm_blocking_space blocking(*m, *n, *m, 1, + false); func[code](*m, *n, a, *lda, b, 1, *ldb, blocking); } else { - internal::gemm_blocking_space blocking(*m, *n, *n, 1, - false); + Eigen::internal::gemm_blocking_space blocking(*m, *n, *n, 1, + false); func[code](*n, *m, a, *lda, b, 1, *ldb, blocking); } @@ -216,89 +242,96 @@ EIGEN_BLAS_FUNC(trmm) const RealScalar *palpha, const RealScalar *pa, const int *lda, RealScalar *pb, const int *ldb) { // std::cerr << "in trmm " << *side << " " << *uplo << " " << *opa << " " << *diag << " " << *m << " " << *n << " " // << *lda << " " << *ldb << " " << *palpha << "\n"; + using Eigen::ColMajor; + using Eigen::DenseIndex; + using Eigen::Dynamic; + using Eigen::Lower; + using Eigen::RowMajor; + using Eigen::UnitDiag; + using Eigen::Upper; typedef void (*functype)(DenseIndex, DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, DenseIndex, const Scalar &, - internal::level3_blocking &); + Eigen::internal::level3_blocking &); static const functype func[32] = { // array index: NOTR | (LEFT << 2) | (UP << 3) | (NUNIT << 4) - (internal::product_triangular_matrix_matrix::run), + (Eigen::internal::product_triangular_matrix_matrix::run), // array index: TR | (LEFT << 2) | (UP << 3) | (NUNIT << 4) - (internal::product_triangular_matrix_matrix::run), + (Eigen::internal::product_triangular_matrix_matrix::run), // array index: ADJ | (LEFT << 2) | (UP << 3) | (NUNIT << 4) - (internal::product_triangular_matrix_matrix::run), + (Eigen::internal::product_triangular_matrix_matrix::run), 0, // array index: NOTR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4) - (internal::product_triangular_matrix_matrix::run), + (Eigen::internal::product_triangular_matrix_matrix::run), // array index: TR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4) - (internal::product_triangular_matrix_matrix::run), + (Eigen::internal::product_triangular_matrix_matrix::run), // array index: ADJ | (RIGHT << 2) | (UP << 3) | (NUNIT << 4) - (internal::product_triangular_matrix_matrix::run), + (Eigen::internal::product_triangular_matrix_matrix::run), 0, // array index: NOTR | (LEFT << 2) | (LO << 3) | (NUNIT << 4) - (internal::product_triangular_matrix_matrix::run), + (Eigen::internal::product_triangular_matrix_matrix::run), // array index: TR | (LEFT << 2) | (LO << 3) | (NUNIT << 4) - (internal::product_triangular_matrix_matrix::run), + (Eigen::internal::product_triangular_matrix_matrix::run), // array index: ADJ | (LEFT << 2) | (LO << 3) | (NUNIT << 4) - (internal::product_triangular_matrix_matrix::run), + (Eigen::internal::product_triangular_matrix_matrix::run), 0, // array index: NOTR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4) - (internal::product_triangular_matrix_matrix::run), + (Eigen::internal::product_triangular_matrix_matrix::run), // array index: TR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4) - (internal::product_triangular_matrix_matrix::run), + (Eigen::internal::product_triangular_matrix_matrix::run), // array index: ADJ | (RIGHT << 2) | (LO << 3) | (NUNIT << 4) - (internal::product_triangular_matrix_matrix::run), + (Eigen::internal::product_triangular_matrix_matrix::run), 0, // array index: NOTR | (LEFT << 2) | (UP << 3) | (UNIT << 4) - (internal::product_triangular_matrix_matrix::run), + (Eigen::internal::product_triangular_matrix_matrix::run), // array index: TR | (LEFT << 2) | (UP << 3) | (UNIT << 4) - (internal::product_triangular_matrix_matrix::run), + (Eigen::internal::product_triangular_matrix_matrix::run), // array index: ADJ | (LEFT << 2) | (UP << 3) | (UNIT << 4) - (internal::product_triangular_matrix_matrix::run), + (Eigen::internal::product_triangular_matrix_matrix::run), 0, // array index: NOTR | (RIGHT << 2) | (UP << 3) | (UNIT << 4) - (internal::product_triangular_matrix_matrix::run), + (Eigen::internal::product_triangular_matrix_matrix::run), // array index: TR | (RIGHT << 2) | (UP << 3) | (UNIT << 4) - (internal::product_triangular_matrix_matrix::run), + (Eigen::internal::product_triangular_matrix_matrix::run), // array index: ADJ | (RIGHT << 2) | (UP << 3) | (UNIT << 4) - (internal::product_triangular_matrix_matrix::run), + (Eigen::internal::product_triangular_matrix_matrix::run), 0, // array index: NOTR | (LEFT << 2) | (LO << 3) | (UNIT << 4) - (internal::product_triangular_matrix_matrix::run), + (Eigen::internal::product_triangular_matrix_matrix::run), // array index: TR | (LEFT << 2) | (LO << 3) | (UNIT << 4) - (internal::product_triangular_matrix_matrix::run), + (Eigen::internal::product_triangular_matrix_matrix::run), // array index: ADJ | (LEFT << 2) | (LO << 3) | (UNIT << 4) - (internal::product_triangular_matrix_matrix::run), + (Eigen::internal::product_triangular_matrix_matrix::run), 0, // array index: NOTR | (RIGHT << 2) | (LO << 3) | (UNIT << 4) - (internal::product_triangular_matrix_matrix::run), + (Eigen::internal::product_triangular_matrix_matrix::run), // array index: TR | (RIGHT << 2) | (LO << 3) | (UNIT << 4) - (internal::product_triangular_matrix_matrix::run), + (Eigen::internal::product_triangular_matrix_matrix::run), // array index: ADJ | (RIGHT << 2) | (LO << 3) | (UNIT << 4) - (internal::product_triangular_matrix_matrix::run), + (Eigen::internal::product_triangular_matrix_matrix::run), 0}; const Scalar *a = reinterpret_cast(pa); @@ -329,16 +362,16 @@ EIGEN_BLAS_FUNC(trmm) if (*m == 0 || *n == 0) return; // FIXME find a way to avoid this copy - Matrix tmp = matrix(b, *m, *n, *ldb); + Eigen::Matrix tmp = matrix(b, *m, *n, *ldb); matrix(b, *m, *n, *ldb).setZero(); if (SIDE(*side) == LEFT) { - internal::gemm_blocking_space blocking(*m, *n, *m, 1, - false); + Eigen::internal::gemm_blocking_space blocking(*m, *n, *m, 1, + false); func[code](*m, *n, *m, a, *lda, tmp.data(), tmp.outerStride(), b, 1, *ldb, alpha, blocking); } else { - internal::gemm_blocking_space blocking(*m, *n, *n, 1, - false); + Eigen::internal::gemm_blocking_space blocking(*m, *n, *n, 1, + false); func[code](*m, *n, *n, tmp.data(), tmp.outerStride(), a, *lda, b, 1, *ldb, alpha, blocking); } } @@ -383,9 +416,15 @@ EIGEN_BLAS_FUNC(symm) if (*m == 0 || *n == 0) return; int size = (SIDE(*side) == LEFT) ? (*m) : (*n); + using Eigen::ColMajor; + using Eigen::DenseIndex; + using Eigen::Dynamic; + using Eigen::Lower; + using Eigen::RowMajor; + using Eigen::Upper; #if ISCOMPLEX // FIXME add support for symmetric complex matrix - Matrix matA(size, size); + Eigen::Matrix matA(size, size); if (UPLO(*uplo) == UP) { matA.triangularView() = matrix(a, size, size, *lda); matA.triangularView() = matrix(a, size, size, *lda).transpose(); @@ -398,24 +437,29 @@ EIGEN_BLAS_FUNC(symm) else if (SIDE(*side) == RIGHT) matrix(c, *m, *n, *ldc) += alpha * matrix(b, *m, *n, *ldb) * matA; #else - internal::gemm_blocking_space blocking(*m, *n, size, 1, false); + Eigen::internal::gemm_blocking_space blocking(*m, *n, size, 1, + false); if (SIDE(*side) == LEFT) if (UPLO(*uplo) == UP) - internal::product_selfadjoint_matrix::run(*m, *n, a, *lda, b, *ldb, c, 1, *ldc, alpha, blocking); + Eigen::internal::product_selfadjoint_matrix::run(*m, *n, a, *lda, b, *ldb, c, 1, *ldc, alpha, + blocking); else if (UPLO(*uplo) == LO) - internal::product_selfadjoint_matrix::run(*m, *n, a, *lda, b, *ldb, c, 1, *ldc, alpha, blocking); + Eigen::internal::product_selfadjoint_matrix::run(*m, *n, a, *lda, b, *ldb, c, 1, *ldc, alpha, + blocking); else return; else if (SIDE(*side) == RIGHT) if (UPLO(*uplo) == UP) - internal::product_selfadjoint_matrix::run(*m, *n, b, *ldb, a, *lda, c, 1, *ldc, alpha, blocking); + Eigen::internal::product_selfadjoint_matrix::run(*m, *n, b, *ldb, a, *lda, c, 1, *ldc, alpha, + blocking); else if (UPLO(*uplo) == LO) - internal::product_selfadjoint_matrix::run(*m, *n, b, *ldb, a, *lda, c, 1, *ldc, alpha, blocking); + Eigen::internal::product_selfadjoint_matrix::run(*m, *n, b, *ldb, a, *lda, c, 1, *ldc, alpha, + blocking); else return; else @@ -430,29 +474,35 @@ EIGEN_BLAS_FUNC(syrk) const int *lda, const RealScalar *pbeta, RealScalar *pc, const int *ldc) { // std::cerr << "in syrk " << *uplo << " " << *op << " " << *n << " " << *k << " " << *palpha << " " << *lda << " " // << *pbeta << " " << *ldc << "\n"; + using Eigen::ColMajor; + using Eigen::DenseIndex; + using Eigen::Dynamic; + using Eigen::Lower; + using Eigen::RowMajor; + using Eigen::Upper; #if !ISCOMPLEX typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, - DenseIndex, DenseIndex, const Scalar &, internal::level3_blocking &); + DenseIndex, DenseIndex, const Scalar &, Eigen::internal::level3_blocking &); static const functype func[8] = { // array index: NOTR | (UP << 2) - (internal::general_matrix_matrix_triangular_product::run), + (Eigen::internal::general_matrix_matrix_triangular_product::run), // array index: TR | (UP << 2) - (internal::general_matrix_matrix_triangular_product::run), + (Eigen::internal::general_matrix_matrix_triangular_product::run), // array index: ADJ | (UP << 2) - (internal::general_matrix_matrix_triangular_product::run), + (Eigen::internal::general_matrix_matrix_triangular_product::run), 0, // array index: NOTR | (LO << 2) - (internal::general_matrix_matrix_triangular_product::run), + (Eigen::internal::general_matrix_matrix_triangular_product::run), // array index: TR | (LO << 2) - (internal::general_matrix_matrix_triangular_product::run), + (Eigen::internal::general_matrix_matrix_triangular_product::run), // array index: ADJ | (LO << 2) - (internal::general_matrix_matrix_triangular_product::run), + (Eigen::internal::general_matrix_matrix_triangular_product::run), 0}; #endif @@ -508,7 +558,8 @@ EIGEN_BLAS_FUNC(syrk) alpha * matrix(a, *k, *n, *lda).transpose() * matrix(a, *k, *n, *lda); } #else - internal::gemm_blocking_space blocking(*n, *n, *k, 1, false); + Eigen::internal::gemm_blocking_space blocking(*n, *n, *k, 1, + false); int code = OP(*op) | (UPLO(*uplo) << 2); func[code](*n, *k, a, *lda, a, *lda, c, 1, *ldc, alpha, blocking); @@ -546,6 +597,8 @@ EIGEN_BLAS_FUNC(syr2k) info = 12; if (info) return xerbla_(SCALAR_SUFFIX_UP "SYR2K", &info); + using Eigen::Lower; + using Eigen::Upper; if (beta != Scalar(1)) { if (UPLO(*uplo) == UP) if (beta == Scalar(0)) @@ -621,16 +674,25 @@ EIGEN_BLAS_FUNC(hemm) if (*m == 0 || *n == 0) return; + using Eigen::ColMajor; + using Eigen::DenseIndex; + using Eigen::Dynamic; + using Eigen::RowMajor; + using Eigen::Upper; + int size = (SIDE(*side) == LEFT) ? (*m) : (*n); - internal::gemm_blocking_space blocking(*m, *n, size, 1, false); + Eigen::internal::gemm_blocking_space blocking(*m, *n, size, 1, + false); if (SIDE(*side) == LEFT) { if (UPLO(*uplo) == UP) - internal::product_selfadjoint_matrix::run(*m, *n, a, *lda, b, *ldb, c, 1, *ldc, alpha, blocking); + Eigen::internal::product_selfadjoint_matrix::run(*m, *n, a, *lda, b, *ldb, c, 1, *ldc, alpha, + blocking); else if (UPLO(*uplo) == LO) - internal::product_selfadjoint_matrix::run(*m, *n, a, *lda, b, *ldb, c, 1, *ldc, alpha, blocking); + Eigen::internal::product_selfadjoint_matrix::run(*m, *n, a, *lda, b, *ldb, c, 1, *ldc, alpha, + blocking); else return; } else if (SIDE(*side) == RIGHT) { @@ -642,8 +704,9 @@ EIGEN_BLAS_FUNC(hemm) RowMajor,true,Conj, ColMajor, 1> ::run(*m, *n, b, *ldb, a, *lda, c, 1, *ldc, alpha, blocking);*/ else if (UPLO(*uplo) == LO) - internal::product_selfadjoint_matrix::run(*m, *n, b, *ldb, a, *lda, c, 1, *ldc, alpha, blocking); + Eigen::internal::product_selfadjoint_matrix::run(*m, *n, b, *ldb, a, *lda, c, 1, *ldc, alpha, + blocking); else return; } else { @@ -658,25 +721,32 @@ EIGEN_BLAS_FUNC(herk) const int *lda, const RealScalar *pbeta, RealScalar *pc, const int *ldc) { // std::cerr << "in herk " << *uplo << " " << *op << " " << *n << " " << *k << " " << *palpha << " " << *lda << " " // << *pbeta << " " << *ldc << "\n"; - + using Eigen::ColMajor; + using Eigen::DenseIndex; + using Eigen::Dynamic; + using Eigen::Lower; + using Eigen::RowMajor; + using Eigen::StrictlyLower; + using Eigen::StrictlyUpper; + using Eigen::Upper; typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, DenseIndex, const Scalar &, Eigen::internal::level3_blocking &); static const functype func[8] = { // array index: NOTR | (UP << 2) - (internal::general_matrix_matrix_triangular_product::run), + (Eigen::internal::general_matrix_matrix_triangular_product::run), 0, // array index: ADJ | (UP << 2) - (internal::general_matrix_matrix_triangular_product::run), + (Eigen::internal::general_matrix_matrix_triangular_product::run), 0, // array index: NOTR | (LO << 2) - (internal::general_matrix_matrix_triangular_product::run), + (Eigen::internal::general_matrix_matrix_triangular_product::run), 0, // array index: ADJ | (LO << 2) - (internal::general_matrix_matrix_triangular_product::run), + (Eigen::internal::general_matrix_matrix_triangular_product::run), 0}; const Scalar *a = reinterpret_cast(pa); @@ -722,7 +792,8 @@ EIGEN_BLAS_FUNC(herk) } if (*k > 0 && alpha != RealScalar(0)) { - internal::gemm_blocking_space blocking(*n, *n, *k, 1, false); + Eigen::internal::gemm_blocking_space blocking(*n, *n, *k, 1, + false); func[code](*n, *k, a, *lda, a, *lda, c, 1, *ldc, alpha, blocking); matrix(c, *n, *n, *ldc).diagonal().imag().setZero(); } @@ -759,6 +830,10 @@ EIGEN_BLAS_FUNC(her2k) info = 12; if (info) return xerbla_(SCALAR_SUFFIX_UP "HER2K", &info); + using Eigen::Lower; + using Eigen::StrictlyLower; + using Eigen::StrictlyUpper; + using Eigen::Upper; if (beta != RealScalar(1)) { if (UPLO(*uplo) == UP) if (beta == Scalar(0)) @@ -783,20 +858,20 @@ EIGEN_BLAS_FUNC(her2k) if (UPLO(*uplo) == UP) { matrix(c, *n, *n, *ldc).triangularView() += alpha * matrix(a, *n, *k, *lda) * matrix(b, *n, *k, *ldb).adjoint() + - numext::conj(alpha) * matrix(b, *n, *k, *ldb) * matrix(a, *n, *k, *lda).adjoint(); + Eigen::numext::conj(alpha) * matrix(b, *n, *k, *ldb) * matrix(a, *n, *k, *lda).adjoint(); } else if (UPLO(*uplo) == LO) matrix(c, *n, *n, *ldc).triangularView() += alpha * matrix(a, *n, *k, *lda) * matrix(b, *n, *k, *ldb).adjoint() + - numext::conj(alpha) * matrix(b, *n, *k, *ldb) * matrix(a, *n, *k, *lda).adjoint(); + Eigen::numext::conj(alpha) * matrix(b, *n, *k, *ldb) * matrix(a, *n, *k, *lda).adjoint(); } else if (OP(*op) == ADJ) { if (UPLO(*uplo) == UP) matrix(c, *n, *n, *ldc).triangularView() += alpha * matrix(a, *k, *n, *lda).adjoint() * matrix(b, *k, *n, *ldb) + - numext::conj(alpha) * matrix(b, *k, *n, *ldb).adjoint() * matrix(a, *k, *n, *lda); + Eigen::numext::conj(alpha) * matrix(b, *k, *n, *ldb).adjoint() * matrix(a, *k, *n, *lda); else if (UPLO(*uplo) == LO) matrix(c, *n, *n, *ldc).triangularView() += alpha * matrix(a, *k, *n, *lda).adjoint() * matrix(b, *k, *n, *ldb) + - numext::conj(alpha) * matrix(b, *k, *n, *ldb).adjoint() * matrix(a, *k, *n, *lda); + Eigen::numext::conj(alpha) * matrix(b, *k, *n, *ldb).adjoint() * matrix(a, *k, *n, *lda); } } diff --git a/lapack/cholesky.inc b/lapack/cholesky.inc index dea5bf635..a93a51152 100644 --- a/lapack/cholesky.inc +++ b/lapack/cholesky.inc @@ -28,9 +28,9 @@ EIGEN_LAPACK_FUNC(potrf)(char *uplo, int *n, RealScalar *pa, int *lda, int *info MatrixType A(a, *n, *n, *lda); int ret; if (UPLO(*uplo) == UP) - ret = int(internal::llt_inplace::blocked(A)); + ret = int(Eigen::internal::llt_inplace::blocked(A)); else - ret = int(internal::llt_inplace::blocked(A)); + ret = int(Eigen::internal::llt_inplace::blocked(A)); if (ret >= 0) *info = ret + 1; } @@ -61,10 +61,10 @@ EIGEN_LAPACK_FUNC(potrs)(char *uplo, int *n, int *nrhs, RealScalar *pa, int *lda MatrixType B(b, *n, *nrhs, *ldb); if (UPLO(*uplo) == UP) { - A.triangularView().adjoint().solveInPlace(B); - A.triangularView().solveInPlace(B); + A.triangularView().adjoint().solveInPlace(B); + A.triangularView().solveInPlace(B); } else { - A.triangularView().solveInPlace(B); - A.triangularView().adjoint().solveInPlace(B); + A.triangularView().solveInPlace(B); + A.triangularView().adjoint().solveInPlace(B); } } diff --git a/lapack/eigenvalues.inc b/lapack/eigenvalues.inc index 6f168de1a..211a7ff4e 100644 --- a/lapack/eigenvalues.inc +++ b/lapack/eigenvalues.inc @@ -47,9 +47,10 @@ EIGEN_LAPACK_FUNC(syev) mat = matrix(a, *n, *n, *lda); bool computeVectors = *jobz == 'V' || *jobz == 'v'; - SelfAdjointEigenSolver eig(mat, computeVectors ? ComputeEigenvectors : EigenvaluesOnly); + Eigen::SelfAdjointEigenSolver eig( + mat, computeVectors ? Eigen::ComputeEigenvectors : Eigen::EigenvaluesOnly); - if (eig.info() == NoConvergence) { + if (eig.info() == Eigen::NoConvergence) { make_vector(w, *n).setZero(); if (computeVectors) matrix(a, *n, *n, *lda).setIdentity(); //*info = 1; diff --git a/lapack/lu.inc b/lapack/lu.inc index d30c8cec6..2ddaf957a 100644 --- a/lapack/lu.inc +++ b/lapack/lu.inc @@ -62,6 +62,8 @@ EIGEN_LAPACK_FUNC(getrs) MatrixType lu(a, *n, *n, *lda); MatrixType B(b, *n, *nrhs, *ldb); + using Eigen::UnitLower; + using Eigen::Upper; for (int i = 0; i < *n; ++i) ipiv[i]--; if (OP(*trans) == NOTR) { B = PivotsType(ipiv, *n) * B; diff --git a/lapack/svd.inc b/lapack/svd.inc index 8e453109f..262c5c696 100644 --- a/lapack/svd.inc +++ b/lapack/svd.inc @@ -56,12 +56,12 @@ EIGEN_LAPACK_FUNC(gesdd) PlainMatrixType mat(*m, *n); mat = matrix(a, *m, *n, *lda); - int option = *jobz == 'A' ? ComputeFullU | ComputeFullV - : *jobz == 'S' ? ComputeThinU | ComputeThinV - : *jobz == 'O' ? ComputeThinU | ComputeThinV + int option = *jobz == 'A' ? Eigen::ComputeFullU | Eigen::ComputeFullV + : *jobz == 'S' ? Eigen::ComputeThinU | Eigen::ComputeThinV + : *jobz == 'O' ? Eigen::ComputeThinU | Eigen::ComputeThinV : 0; - BDCSVD svd(mat, option); + Eigen::BDCSVD svd(mat, option); make_vector(s, diag_size) = svd.singularValues().head(diag_size); @@ -119,14 +119,14 @@ EIGEN_LAPACK_FUNC(gesvd) PlainMatrixType mat(*m, *n); mat = matrix(a, *m, *n, *lda); - int option = (*jobu == 'A' ? ComputeFullU - : *jobu == 'S' || *jobu == 'O' ? ComputeThinU + int option = (*jobu == 'A' ? Eigen::ComputeFullU + : *jobu == 'S' || *jobu == 'O' ? Eigen::ComputeThinU : 0) | - (*jobv == 'A' ? ComputeFullV - : *jobv == 'S' || *jobv == 'O' ? ComputeThinV + (*jobv == 'A' ? Eigen::ComputeFullV + : *jobv == 'S' || *jobv == 'O' ? Eigen::ComputeThinV : 0); - JacobiSVD svd(mat, option); + Eigen::JacobiSVD svd(mat, option); make_vector(s, diag_size) = svd.singularValues().head(diag_size); {