| // 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" | |
| EIGEN_BLAS_FUNC(gemm) | |
| (const char *opa, const char *opb, const EIGEN_BLAS_INT *m, const EIGEN_BLAS_INT *n, const EIGEN_BLAS_INT *k, | |
| const RealScalar *palpha, const RealScalar *pa, const EIGEN_BLAS_INT *lda, const RealScalar *pb, | |
| const EIGEN_BLAS_INT *ldb, const RealScalar *pbeta, RealScalar *pc, const EIGEN_BLAS_INT *ldc) { | |
| 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<Scalar, Scalar> &, | |
| Eigen::internal::GemmParallelInfo<DenseIndex> *); | |
| static const functype func[12] = { | |
| // array index: NOTR | (NOTR << 2) | |
| (Eigen::internal::general_matrix_matrix_product<DenseIndex, Scalar, ColMajor, false, Scalar, ColMajor, false, | |
| ColMajor, 1>::run), | |
| // array index: TR | (NOTR << 2) | |
| (Eigen::internal::general_matrix_matrix_product<DenseIndex, Scalar, RowMajor, false, Scalar, ColMajor, false, | |
| ColMajor, 1>::run), | |
| // array index: ADJ | (NOTR << 2) | |
| (Eigen::internal::general_matrix_matrix_product<DenseIndex, Scalar, RowMajor, Conj, Scalar, ColMajor, false, | |
| ColMajor, 1>::run), | |
| 0, | |
| // array index: NOTR | (TR << 2) | |
| (Eigen::internal::general_matrix_matrix_product<DenseIndex, Scalar, ColMajor, false, Scalar, RowMajor, false, | |
| ColMajor, 1>::run), | |
| // array index: TR | (TR << 2) | |
| (Eigen::internal::general_matrix_matrix_product<DenseIndex, Scalar, RowMajor, false, Scalar, RowMajor, false, | |
| ColMajor, 1>::run), | |
| // array index: ADJ | (TR << 2) | |
| (Eigen::internal::general_matrix_matrix_product<DenseIndex, Scalar, RowMajor, Conj, Scalar, RowMajor, false, | |
| ColMajor, 1>::run), | |
| 0, | |
| // array index: NOTR | (ADJ << 2) | |
| (Eigen::internal::general_matrix_matrix_product<DenseIndex, Scalar, ColMajor, false, Scalar, RowMajor, Conj, | |
| ColMajor, 1>::run), | |
| // array index: TR | (ADJ << 2) | |
| (Eigen::internal::general_matrix_matrix_product<DenseIndex, Scalar, RowMajor, false, Scalar, RowMajor, Conj, | |
| ColMajor, 1>::run), | |
| // array index: ADJ | (ADJ << 2) | |
| (Eigen::internal::general_matrix_matrix_product<DenseIndex, Scalar, RowMajor, Conj, Scalar, RowMajor, Conj, | |
| ColMajor, 1>::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); | |
| EIGEN_BLAS_INT info = 0; | |
| if (OP(*opa) == INVALID) | |
| info = 1; | |
| else if (OP(*opb) == INVALID) | |
| info = 2; | |
| else if (*m < 0) | |
| info = 3; | |
| else if (*n < 0) | |
| info = 4; | |
| else if (*k < 0) | |
| info = 5; | |
| else if (*lda < std::max<EIGEN_BLAS_INT>(1, (OP(*opa) == NOTR) ? *m : *k)) | |
| info = 8; | |
| else if (*ldb < std::max<EIGEN_BLAS_INT>(1, (OP(*opb) == NOTR) ? *k : *n)) | |
| info = 10; | |
| else if (*ldc < std::max<EIGEN_BLAS_INT>(1, *m)) | |
| info = 13; | |
| if (info) return xerbla_(SCALAR_SUFFIX_UP "GEMM ", &info, kBlasNameLength); | |
| if (*m == 0 || *n == 0) return; | |
| if (beta != Scalar(1)) { | |
| if (beta == Scalar(0)) | |
| matrix(c, *m, *n, *ldc).setZero(); | |
| else | |
| matrix(c, *m, *n, *ldc) *= beta; | |
| } | |
| if (*k == 0) return; | |
| Eigen::internal::gemm_blocking_space<ColMajor, Scalar, Scalar, Dynamic, Dynamic, Dynamic> 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); | |
| } | |
| EIGEN_BLAS_FUNC(gemmtr) | |
| (const char *uplo, const char *opa, const char *opb, const EIGEN_BLAS_INT *n, const EIGEN_BLAS_INT *k, | |
| const RealScalar *palpha, const RealScalar *pa, const EIGEN_BLAS_INT *lda, const RealScalar *pb, | |
| const EIGEN_BLAS_INT *ldb, const RealScalar *pbeta, RealScalar *pc, const EIGEN_BLAS_INT *ldc) { | |
| using Eigen::ColMajor; | |
| using Eigen::DenseIndex; | |
| using Eigen::Dynamic; | |
| using Eigen::Lower; | |
| using Eigen::RowMajor; | |
| using Eigen::Upper; | |
| typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, | |
| DenseIndex, DenseIndex, const Scalar &, Eigen::internal::level3_blocking<Scalar, Scalar> &); | |
| static const functype func[24] = { | |
| // Upper-triangular result. | |
| // array index: NOTR | (NOTR << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, ColMajor, false, Scalar, ColMajor, | |
| false, ColMajor, 1, Upper>::run), | |
| // array index: TR | (NOTR << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, RowMajor, false, Scalar, ColMajor, | |
| false, ColMajor, 1, Upper>::run), | |
| // array index: ADJ | (NOTR << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, RowMajor, Conj, Scalar, ColMajor, | |
| false, ColMajor, 1, Upper>::run), | |
| 0, | |
| // array index: NOTR | (TR << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, ColMajor, false, Scalar, RowMajor, | |
| false, ColMajor, 1, Upper>::run), | |
| // array index: TR | (TR << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, RowMajor, false, Scalar, RowMajor, | |
| false, ColMajor, 1, Upper>::run), | |
| // array index: ADJ | (TR << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, RowMajor, Conj, Scalar, RowMajor, | |
| false, ColMajor, 1, Upper>::run), | |
| 0, | |
| // array index: NOTR | (ADJ << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, ColMajor, false, Scalar, RowMajor, | |
| Conj, ColMajor, 1, Upper>::run), | |
| // array index: TR | (ADJ << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, RowMajor, false, Scalar, RowMajor, | |
| Conj, ColMajor, 1, Upper>::run), | |
| // array index: ADJ | (ADJ << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, RowMajor, Conj, Scalar, RowMajor, | |
| Conj, ColMajor, 1, Upper>::run), | |
| 0, | |
| // Lower-triangular result. | |
| // array index: NOTR | (NOTR << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, ColMajor, false, Scalar, ColMajor, | |
| false, ColMajor, 1, Lower>::run), | |
| // array index: TR | (NOTR << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, RowMajor, false, Scalar, ColMajor, | |
| false, ColMajor, 1, Lower>::run), | |
| // array index: ADJ | (NOTR << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, RowMajor, Conj, Scalar, ColMajor, | |
| false, ColMajor, 1, Lower>::run), | |
| 0, | |
| // array index: NOTR | (TR << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, ColMajor, false, Scalar, RowMajor, | |
| false, ColMajor, 1, Lower>::run), | |
| // array index: TR | (TR << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, RowMajor, false, Scalar, RowMajor, | |
| false, ColMajor, 1, Lower>::run), | |
| // array index: ADJ | (TR << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, RowMajor, Conj, Scalar, RowMajor, | |
| false, ColMajor, 1, Lower>::run), | |
| 0, | |
| // array index: NOTR | (ADJ << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, ColMajor, false, Scalar, RowMajor, | |
| Conj, ColMajor, 1, Lower>::run), | |
| // array index: TR | (ADJ << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, RowMajor, false, Scalar, RowMajor, | |
| Conj, ColMajor, 1, Lower>::run), | |
| // array index: ADJ | (ADJ << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, RowMajor, Conj, Scalar, RowMajor, | |
| Conj, ColMajor, 1, Lower>::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); | |
| EIGEN_BLAS_INT info = 0; | |
| if (UPLO(*uplo) == INVALID) | |
| info = 1; | |
| else if (OP(*opa) == INVALID) | |
| info = 2; | |
| else if (OP(*opb) == INVALID) | |
| info = 3; | |
| else if (*n < 0) | |
| info = 4; | |
| else if (*k < 0) | |
| info = 5; | |
| else if (*lda < std::max<EIGEN_BLAS_INT>(1, (OP(*opa) == NOTR) ? *n : *k)) | |
| info = 8; | |
| else if (*ldb < std::max<EIGEN_BLAS_INT>(1, (OP(*opb) == NOTR) ? *k : *n)) | |
| info = 10; | |
| else if (*ldc < std::max<EIGEN_BLAS_INT>(1, *n)) | |
| info = 13; | |
| if (info) return xerbla_(SCALAR_SUFFIX_UP "GEMMTR", &info, kBlasNameLength + 1); | |
| if (*n == 0) return; | |
| const int upper = (UPLO(*uplo) == UP); | |
| if (beta != Scalar(1)) { | |
| if (beta == Scalar(0)) { | |
| if (upper) { | |
| matrix(c, *n, *n, *ldc).triangularView<Eigen::Upper>().setZero(); | |
| } else { | |
| matrix(c, *n, *n, *ldc).triangularView<Eigen::Lower>().setZero(); | |
| } | |
| } else { | |
| if (upper) { | |
| matrix(c, *n, *n, *ldc).triangularView<Eigen::Upper>() *= beta; | |
| } else { | |
| matrix(c, *n, *n, *ldc).triangularView<Eigen::Lower>() *= beta; | |
| } | |
| } | |
| } | |
| if (*k == 0) return; | |
| Eigen::internal::gemm_blocking_space<ColMajor, Scalar, Scalar, Dynamic, Dynamic, Dynamic> blocking(*n, *n, *k, 1, | |
| false); | |
| int code = OP(*opa) | ((OP(*opb) << 2) + (UPLO(*uplo) * 12)); | |
| func[code](*n, *k, a, *lda, b, *ldb, c, 1, *ldc, alpha, blocking); | |
| } | |
| EIGEN_BLAS_FUNC(trsm) | |
| (const char *side, const char *uplo, const char *opa, const char *diag, const EIGEN_BLAS_INT *m, | |
| const EIGEN_BLAS_INT *n, const RealScalar *palpha, const RealScalar *pa, const EIGEN_BLAS_INT *lda, RealScalar *pb, | |
| const EIGEN_BLAS_INT *ldb) { | |
| 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<Scalar, Scalar> &); | |
| static const functype func[32] = { | |
| // array index: NOTR | (LEFT << 2) | (UP << 3) | (NUNIT << 4) | |
| (Eigen::internal::triangular_solve_matrix<Scalar, DenseIndex, OnTheLeft, Upper | 0, false, ColMajor, ColMajor, | |
| 1>::run), | |
| // array index: TR | (LEFT << 2) | (UP << 3) | (NUNIT << 4) | |
| (Eigen::internal::triangular_solve_matrix<Scalar, DenseIndex, OnTheLeft, Lower | 0, false, RowMajor, ColMajor, | |
| 1>::run), | |
| // array index: ADJ | (LEFT << 2) | (UP << 3) | (NUNIT << 4) | |
| (Eigen::internal::triangular_solve_matrix<Scalar, DenseIndex, OnTheLeft, Lower | 0, Conj, RowMajor, ColMajor, | |
| 1>::run), | |
| 0, | |
| // array index: NOTR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4) | |
| (Eigen::internal::triangular_solve_matrix<Scalar, DenseIndex, OnTheRight, Upper | 0, false, ColMajor, ColMajor, | |
| 1>::run), | |
| // array index: TR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4) | |
| (Eigen::internal::triangular_solve_matrix<Scalar, DenseIndex, OnTheRight, Lower | 0, false, RowMajor, ColMajor, | |
| 1>::run), | |
| // array index: ADJ | (RIGHT << 2) | (UP << 3) | (NUNIT << 4) | |
| (Eigen::internal::triangular_solve_matrix<Scalar, DenseIndex, OnTheRight, Lower | 0, Conj, RowMajor, ColMajor, | |
| 1>::run), | |
| 0, | |
| // array index: NOTR | (LEFT << 2) | (LO << 3) | (NUNIT << 4) | |
| (Eigen::internal::triangular_solve_matrix<Scalar, DenseIndex, OnTheLeft, Lower | 0, false, ColMajor, ColMajor, | |
| 1>::run), | |
| // array index: TR | (LEFT << 2) | (LO << 3) | (NUNIT << 4) | |
| (Eigen::internal::triangular_solve_matrix<Scalar, DenseIndex, OnTheLeft, Upper | 0, false, RowMajor, ColMajor, | |
| 1>::run), | |
| // array index: ADJ | (LEFT << 2) | (LO << 3) | (NUNIT << 4) | |
| (Eigen::internal::triangular_solve_matrix<Scalar, DenseIndex, OnTheLeft, Upper | 0, Conj, RowMajor, ColMajor, | |
| 1>::run), | |
| 0, | |
| // array index: NOTR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4) | |
| (Eigen::internal::triangular_solve_matrix<Scalar, DenseIndex, OnTheRight, Lower | 0, false, ColMajor, ColMajor, | |
| 1>::run), | |
| // array index: TR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4) | |
| (Eigen::internal::triangular_solve_matrix<Scalar, DenseIndex, OnTheRight, Upper | 0, false, RowMajor, ColMajor, | |
| 1>::run), | |
| // array index: ADJ | (RIGHT << 2) | (LO << 3) | (NUNIT << 4) | |
| (Eigen::internal::triangular_solve_matrix<Scalar, DenseIndex, OnTheRight, Upper | 0, Conj, RowMajor, ColMajor, | |
| 1>::run), | |
| 0, | |
| // array index: NOTR | (LEFT << 2) | (UP << 3) | (UNIT << 4) | |
| (Eigen::internal::triangular_solve_matrix<Scalar, DenseIndex, OnTheLeft, Upper | UnitDiag, false, ColMajor, | |
| ColMajor, 1>::run), | |
| // array index: TR | (LEFT << 2) | (UP << 3) | (UNIT << 4) | |
| (Eigen::internal::triangular_solve_matrix<Scalar, DenseIndex, OnTheLeft, Lower | UnitDiag, false, RowMajor, | |
| ColMajor, 1>::run), | |
| // array index: ADJ | (LEFT << 2) | (UP << 3) | (UNIT << 4) | |
| (Eigen::internal::triangular_solve_matrix<Scalar, DenseIndex, OnTheLeft, Lower | UnitDiag, Conj, RowMajor, | |
| ColMajor, 1>::run), | |
| 0, | |
| // array index: NOTR | (RIGHT << 2) | (UP << 3) | (UNIT << 4) | |
| (Eigen::internal::triangular_solve_matrix<Scalar, DenseIndex, OnTheRight, Upper | UnitDiag, false, ColMajor, | |
| ColMajor, 1>::run), | |
| // array index: TR | (RIGHT << 2) | (UP << 3) | (UNIT << 4) | |
| (Eigen::internal::triangular_solve_matrix<Scalar, DenseIndex, OnTheRight, Lower | UnitDiag, false, RowMajor, | |
| ColMajor, 1>::run), | |
| // array index: ADJ | (RIGHT << 2) | (UP << 3) | (UNIT << 4) | |
| (Eigen::internal::triangular_solve_matrix<Scalar, DenseIndex, OnTheRight, Lower | UnitDiag, Conj, RowMajor, | |
| ColMajor, 1>::run), | |
| 0, | |
| // array index: NOTR | (LEFT << 2) | (LO << 3) | (UNIT << 4) | |
| (Eigen::internal::triangular_solve_matrix<Scalar, DenseIndex, OnTheLeft, Lower | UnitDiag, false, ColMajor, | |
| ColMajor, 1>::run), | |
| // array index: TR | (LEFT << 2) | (LO << 3) | (UNIT << 4) | |
| (Eigen::internal::triangular_solve_matrix<Scalar, DenseIndex, OnTheLeft, Upper | UnitDiag, false, RowMajor, | |
| ColMajor, 1>::run), | |
| // array index: ADJ | (LEFT << 2) | (LO << 3) | (UNIT << 4) | |
| (Eigen::internal::triangular_solve_matrix<Scalar, DenseIndex, OnTheLeft, Upper | UnitDiag, Conj, RowMajor, | |
| ColMajor, 1>::run), | |
| 0, | |
| // array index: NOTR | (RIGHT << 2) | (LO << 3) | (UNIT << 4) | |
| (Eigen::internal::triangular_solve_matrix<Scalar, DenseIndex, OnTheRight, Lower | UnitDiag, false, ColMajor, | |
| ColMajor, 1>::run), | |
| // array index: TR | (RIGHT << 2) | (LO << 3) | (UNIT << 4) | |
| (Eigen::internal::triangular_solve_matrix<Scalar, DenseIndex, OnTheRight, Upper | UnitDiag, false, RowMajor, | |
| ColMajor, 1>::run), | |
| // array index: ADJ | (RIGHT << 2) | (LO << 3) | (UNIT << 4) | |
| (Eigen::internal::triangular_solve_matrix<Scalar, DenseIndex, OnTheRight, Upper | UnitDiag, Conj, RowMajor, | |
| ColMajor, 1>::run), | |
| 0}; | |
| const Scalar *a = reinterpret_cast<const Scalar *>(pa); | |
| Scalar *b = reinterpret_cast<Scalar *>(pb); | |
| Scalar alpha = *reinterpret_cast<const Scalar *>(palpha); | |
| EIGEN_BLAS_INT info = 0; | |
| if (SIDE(*side) == INVALID) | |
| info = 1; | |
| else if (UPLO(*uplo) == INVALID) | |
| info = 2; | |
| else if (OP(*opa) == INVALID) | |
| info = 3; | |
| else if (DIAG(*diag) == INVALID) | |
| info = 4; | |
| else if (*m < 0) | |
| info = 5; | |
| else if (*n < 0) | |
| info = 6; | |
| else if (*lda < std::max<EIGEN_BLAS_INT>(1, (SIDE(*side) == LEFT) ? *m : *n)) | |
| info = 9; | |
| else if (*ldb < std::max<EIGEN_BLAS_INT>(1, *m)) | |
| info = 11; | |
| if (info) return xerbla_(SCALAR_SUFFIX_UP "TRSM ", &info, kBlasNameLength); | |
| if (*m == 0 || *n == 0) return; | |
| // When alpha is zero, A is not referenced and B need not be set on entry. | |
| if (alpha == Scalar(0)) { | |
| matrix(b, *m, *n, *ldb).setZero(); | |
| return; | |
| } | |
| int code = OP(*opa) | (SIDE(*side) << 2) | (UPLO(*uplo) << 3) | (DIAG(*diag) << 4); | |
| if (SIDE(*side) == LEFT) { | |
| Eigen::internal::gemm_blocking_space<ColMajor, Scalar, Scalar, Dynamic, Dynamic, Dynamic, 4> blocking(*m, *n, *m, 1, | |
| false); | |
| func[code](*m, *n, a, *lda, b, 1, *ldb, blocking); | |
| } else { | |
| Eigen::internal::gemm_blocking_space<ColMajor, Scalar, Scalar, Dynamic, Dynamic, Dynamic, 4> blocking(*m, *n, *n, 1, | |
| false); | |
| func[code](*n, *m, a, *lda, b, 1, *ldb, blocking); | |
| } | |
| if (alpha != Scalar(1)) matrix(b, *m, *n, *ldb) *= alpha; | |
| } | |
| // b = alpha*op(a)*b for side = 'L'or'l' | |
| // b = alpha*b*op(a) for side = 'R'or'r' | |
| EIGEN_BLAS_FUNC(trmm) | |
| (const char *side, const char *uplo, const char *opa, const char *diag, const EIGEN_BLAS_INT *m, | |
| const EIGEN_BLAS_INT *n, const RealScalar *palpha, const RealScalar *pa, const EIGEN_BLAS_INT *lda, RealScalar *pb, | |
| const EIGEN_BLAS_INT *ldb) { | |
| 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 &, | |
| Eigen::internal::level3_blocking<Scalar, Scalar> &); | |
| static const functype func[32] = { | |
| // array index: NOTR | (LEFT << 2) | (UP << 3) | (NUNIT << 4) | |
| (Eigen::internal::product_triangular_matrix_matrix<Scalar, DenseIndex, Upper | 0, true, ColMajor, false, ColMajor, | |
| false, ColMajor, 1>::run), | |
| // array index: TR | (LEFT << 2) | (UP << 3) | (NUNIT << 4) | |
| (Eigen::internal::product_triangular_matrix_matrix<Scalar, DenseIndex, Lower | 0, true, RowMajor, false, ColMajor, | |
| false, ColMajor, 1>::run), | |
| // array index: ADJ | (LEFT << 2) | (UP << 3) | (NUNIT << 4) | |
| (Eigen::internal::product_triangular_matrix_matrix<Scalar, DenseIndex, Lower | 0, true, RowMajor, Conj, ColMajor, | |
| false, ColMajor, 1>::run), | |
| 0, | |
| // array index: NOTR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4) | |
| (Eigen::internal::product_triangular_matrix_matrix<Scalar, DenseIndex, Upper | 0, false, ColMajor, false, | |
| ColMajor, false, ColMajor, 1>::run), | |
| // array index: TR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4) | |
| (Eigen::internal::product_triangular_matrix_matrix<Scalar, DenseIndex, Lower | 0, false, ColMajor, false, | |
| RowMajor, false, ColMajor, 1>::run), | |
| // array index: ADJ | (RIGHT << 2) | (UP << 3) | (NUNIT << 4) | |
| (Eigen::internal::product_triangular_matrix_matrix<Scalar, DenseIndex, Lower | 0, false, ColMajor, false, | |
| RowMajor, Conj, ColMajor, 1>::run), | |
| 0, | |
| // array index: NOTR | (LEFT << 2) | (LO << 3) | (NUNIT << 4) | |
| (Eigen::internal::product_triangular_matrix_matrix<Scalar, DenseIndex, Lower | 0, true, ColMajor, false, ColMajor, | |
| false, ColMajor, 1>::run), | |
| // array index: TR | (LEFT << 2) | (LO << 3) | (NUNIT << 4) | |
| (Eigen::internal::product_triangular_matrix_matrix<Scalar, DenseIndex, Upper | 0, true, RowMajor, false, ColMajor, | |
| false, ColMajor, 1>::run), | |
| // array index: ADJ | (LEFT << 2) | (LO << 3) | (NUNIT << 4) | |
| (Eigen::internal::product_triangular_matrix_matrix<Scalar, DenseIndex, Upper | 0, true, RowMajor, Conj, ColMajor, | |
| false, ColMajor, 1>::run), | |
| 0, | |
| // array index: NOTR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4) | |
| (Eigen::internal::product_triangular_matrix_matrix<Scalar, DenseIndex, Lower | 0, false, ColMajor, false, | |
| ColMajor, false, ColMajor, 1>::run), | |
| // array index: TR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4) | |
| (Eigen::internal::product_triangular_matrix_matrix<Scalar, DenseIndex, Upper | 0, false, ColMajor, false, | |
| RowMajor, false, ColMajor, 1>::run), | |
| // array index: ADJ | (RIGHT << 2) | (LO << 3) | (NUNIT << 4) | |
| (Eigen::internal::product_triangular_matrix_matrix<Scalar, DenseIndex, Upper | 0, false, ColMajor, false, | |
| RowMajor, Conj, ColMajor, 1>::run), | |
| 0, | |
| // array index: NOTR | (LEFT << 2) | (UP << 3) | (UNIT << 4) | |
| (Eigen::internal::product_triangular_matrix_matrix<Scalar, DenseIndex, Upper | UnitDiag, true, ColMajor, false, | |
| ColMajor, false, ColMajor, 1>::run), | |
| // array index: TR | (LEFT << 2) | (UP << 3) | (UNIT << 4) | |
| (Eigen::internal::product_triangular_matrix_matrix<Scalar, DenseIndex, Lower | UnitDiag, true, RowMajor, false, | |
| ColMajor, false, ColMajor, 1>::run), | |
| // array index: ADJ | (LEFT << 2) | (UP << 3) | (UNIT << 4) | |
| (Eigen::internal::product_triangular_matrix_matrix<Scalar, DenseIndex, Lower | UnitDiag, true, RowMajor, Conj, | |
| ColMajor, false, ColMajor, 1>::run), | |
| 0, | |
| // array index: NOTR | (RIGHT << 2) | (UP << 3) | (UNIT << 4) | |
| (Eigen::internal::product_triangular_matrix_matrix<Scalar, DenseIndex, Upper | UnitDiag, false, ColMajor, false, | |
| ColMajor, false, ColMajor, 1>::run), | |
| // array index: TR | (RIGHT << 2) | (UP << 3) | (UNIT << 4) | |
| (Eigen::internal::product_triangular_matrix_matrix<Scalar, DenseIndex, Lower | UnitDiag, false, ColMajor, false, | |
| RowMajor, false, ColMajor, 1>::run), | |
| // array index: ADJ | (RIGHT << 2) | (UP << 3) | (UNIT << 4) | |
| (Eigen::internal::product_triangular_matrix_matrix<Scalar, DenseIndex, Lower | UnitDiag, false, ColMajor, false, | |
| RowMajor, Conj, ColMajor, 1>::run), | |
| 0, | |
| // array index: NOTR | (LEFT << 2) | (LO << 3) | (UNIT << 4) | |
| (Eigen::internal::product_triangular_matrix_matrix<Scalar, DenseIndex, Lower | UnitDiag, true, ColMajor, false, | |
| ColMajor, false, ColMajor, 1>::run), | |
| // array index: TR | (LEFT << 2) | (LO << 3) | (UNIT << 4) | |
| (Eigen::internal::product_triangular_matrix_matrix<Scalar, DenseIndex, Upper | UnitDiag, true, RowMajor, false, | |
| ColMajor, false, ColMajor, 1>::run), | |
| // array index: ADJ | (LEFT << 2) | (LO << 3) | (UNIT << 4) | |
| (Eigen::internal::product_triangular_matrix_matrix<Scalar, DenseIndex, Upper | UnitDiag, true, RowMajor, Conj, | |
| ColMajor, false, ColMajor, 1>::run), | |
| 0, | |
| // array index: NOTR | (RIGHT << 2) | (LO << 3) | (UNIT << 4) | |
| (Eigen::internal::product_triangular_matrix_matrix<Scalar, DenseIndex, Lower | UnitDiag, false, ColMajor, false, | |
| ColMajor, false, ColMajor, 1>::run), | |
| // array index: TR | (RIGHT << 2) | (LO << 3) | (UNIT << 4) | |
| (Eigen::internal::product_triangular_matrix_matrix<Scalar, DenseIndex, Upper | UnitDiag, false, ColMajor, false, | |
| RowMajor, false, ColMajor, 1>::run), | |
| // array index: ADJ | (RIGHT << 2) | (LO << 3) | (UNIT << 4) | |
| (Eigen::internal::product_triangular_matrix_matrix<Scalar, DenseIndex, Upper | UnitDiag, false, ColMajor, false, | |
| RowMajor, Conj, ColMajor, 1>::run), | |
| 0}; | |
| const Scalar *a = reinterpret_cast<const Scalar *>(pa); | |
| Scalar *b = reinterpret_cast<Scalar *>(pb); | |
| Scalar alpha = *reinterpret_cast<const Scalar *>(palpha); | |
| EIGEN_BLAS_INT info = 0; | |
| if (SIDE(*side) == INVALID) | |
| info = 1; | |
| else if (UPLO(*uplo) == INVALID) | |
| info = 2; | |
| else if (OP(*opa) == INVALID) | |
| info = 3; | |
| else if (DIAG(*diag) == INVALID) | |
| info = 4; | |
| else if (*m < 0) | |
| info = 5; | |
| else if (*n < 0) | |
| info = 6; | |
| else if (*lda < std::max<EIGEN_BLAS_INT>(1, (SIDE(*side) == LEFT) ? *m : *n)) | |
| info = 9; | |
| else if (*ldb < std::max<EIGEN_BLAS_INT>(1, *m)) | |
| info = 11; | |
| if (info) return xerbla_(SCALAR_SUFFIX_UP "TRMM ", &info, kBlasNameLength); | |
| int code = OP(*opa) | (SIDE(*side) << 2) | (UPLO(*uplo) << 3) | (DIAG(*diag) << 4); | |
| if (*m == 0 || *n == 0) return; | |
| // When alpha is zero, A is not referenced and B need not be set on entry. | |
| if (alpha == Scalar(0)) { | |
| matrix(b, *m, *n, *ldb).setZero(); | |
| return; | |
| } | |
| // FIXME: find a way to avoid this copy | |
| Eigen::Matrix<Scalar, Dynamic, Dynamic, ColMajor> tmp = matrix(b, *m, *n, *ldb); | |
| matrix(b, *m, *n, *ldb).setZero(); | |
| if (SIDE(*side) == LEFT) { | |
| Eigen::internal::gemm_blocking_space<ColMajor, Scalar, Scalar, Dynamic, Dynamic, Dynamic, 4> blocking(*m, *n, *m, 1, | |
| false); | |
| func[code](*m, *n, *m, a, *lda, tmp.data(), tmp.outerStride(), b, 1, *ldb, alpha, blocking); | |
| } else { | |
| Eigen::internal::gemm_blocking_space<ColMajor, Scalar, Scalar, Dynamic, Dynamic, Dynamic, 4> blocking(*m, *n, *n, 1, | |
| false); | |
| func[code](*m, *n, *n, tmp.data(), tmp.outerStride(), a, *lda, b, 1, *ldb, alpha, blocking); | |
| } | |
| } | |
| // c = alpha*a*b + beta*c for side = 'L'or'l' | |
| // c = alpha*b*a + beta*c for side = 'R'or'r | |
| EIGEN_BLAS_FUNC(symm) | |
| (const char *side, const char *uplo, 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 *ldb, | |
| const RealScalar *pbeta, RealScalar *pc, const EIGEN_BLAS_INT *ldc) { | |
| 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); | |
| EIGEN_BLAS_INT info = 0; | |
| if (SIDE(*side) == INVALID) | |
| info = 1; | |
| else if (UPLO(*uplo) == INVALID) | |
| info = 2; | |
| else if (*m < 0) | |
| info = 3; | |
| else if (*n < 0) | |
| info = 4; | |
| else if (*lda < std::max<EIGEN_BLAS_INT>(1, (SIDE(*side) == LEFT) ? *m : *n)) | |
| info = 7; | |
| else if (*ldb < std::max<EIGEN_BLAS_INT>(1, *m)) | |
| info = 9; | |
| else if (*ldc < std::max<EIGEN_BLAS_INT>(1, *m)) | |
| info = 12; | |
| if (info) return xerbla_(SCALAR_SUFFIX_UP "SYMM ", &info, kBlasNameLength); | |
| if (beta != Scalar(1)) { | |
| if (beta == Scalar(0)) | |
| matrix(c, *m, *n, *ldc).setZero(); | |
| else | |
| matrix(c, *m, *n, *ldc) *= beta; | |
| } | |
| if (*m == 0 || *n == 0 || alpha == Scalar(0)) return; | |
| EIGEN_BLAS_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 | |
| Eigen::Matrix<Scalar, Dynamic, Dynamic, ColMajor> matA(size, size); | |
| if (UPLO(*uplo) == UP) { | |
| matA.triangularView<Upper>() = matrix(a, size, size, *lda); | |
| matA.triangularView<Lower>() = matrix(a, size, size, *lda).transpose(); | |
| } else if (UPLO(*uplo) == LO) { | |
| matA.triangularView<Lower>() = matrix(a, size, size, *lda); | |
| matA.triangularView<Upper>() = matrix(a, size, size, *lda).transpose(); | |
| } | |
| if (SIDE(*side) == LEFT) | |
| matrix(c, *m, *n, *ldc) += alpha * matA * matrix(b, *m, *n, *ldb); | |
| else if (SIDE(*side) == RIGHT) | |
| matrix(c, *m, *n, *ldc) += alpha * matrix(b, *m, *n, *ldb) * matA; | |
| #else | |
| Eigen::internal::gemm_blocking_space<ColMajor, Scalar, Scalar, Dynamic, Dynamic, Dynamic> blocking(*m, *n, size, 1, | |
| false); | |
| if (SIDE(*side) == LEFT) | |
| if (UPLO(*uplo) == UP) | |
| Eigen::internal::product_selfadjoint_matrix<Scalar, DenseIndex, RowMajor, true, false, ColMajor, false, false, | |
| ColMajor, 1>::run(*m, *n, a, *lda, b, *ldb, c, 1, *ldc, alpha, | |
| blocking); | |
| else if (UPLO(*uplo) == LO) | |
| Eigen::internal::product_selfadjoint_matrix<Scalar, DenseIndex, ColMajor, true, false, ColMajor, false, false, | |
| ColMajor, 1>::run(*m, *n, a, *lda, b, *ldb, c, 1, *ldc, alpha, | |
| blocking); | |
| else | |
| return; | |
| else if (SIDE(*side) == RIGHT) | |
| if (UPLO(*uplo) == UP) | |
| Eigen::internal::product_selfadjoint_matrix<Scalar, DenseIndex, ColMajor, false, false, RowMajor, true, false, | |
| ColMajor, 1>::run(*m, *n, b, *ldb, a, *lda, c, 1, *ldc, alpha, | |
| blocking); | |
| else if (UPLO(*uplo) == LO) | |
| Eigen::internal::product_selfadjoint_matrix<Scalar, DenseIndex, ColMajor, false, false, ColMajor, true, false, | |
| ColMajor, 1>::run(*m, *n, b, *ldb, a, *lda, c, 1, *ldc, alpha, | |
| blocking); | |
| else | |
| return; | |
| else | |
| return; | |
| #endif | |
| } | |
| // c = alpha*a*a' + beta*c for op = 'N'or'n' | |
| // c = alpha*a'*a + beta*c for op = 'T'or't','C'or'c' | |
| EIGEN_BLAS_FUNC(syrk) | |
| (const char *uplo, const char *op, const EIGEN_BLAS_INT *n, const EIGEN_BLAS_INT *k, const RealScalar *palpha, | |
| const RealScalar *pa, const EIGEN_BLAS_INT *lda, const RealScalar *pbeta, RealScalar *pc, const EIGEN_BLAS_INT *ldc) { | |
| 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 &, Eigen::internal::level3_blocking<Scalar, Scalar> &); | |
| static const functype func[8] = { | |
| // array index: NOTR | (UP << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, ColMajor, false, Scalar, RowMajor, | |
| Conj, ColMajor, 1, Upper>::run), | |
| // array index: TR | (UP << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, RowMajor, false, Scalar, ColMajor, | |
| Conj, ColMajor, 1, Upper>::run), | |
| // array index: ADJ | (UP << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, RowMajor, Conj, Scalar, ColMajor, | |
| false, ColMajor, 1, Upper>::run), | |
| 0, | |
| // array index: NOTR | (LO << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, ColMajor, false, Scalar, RowMajor, | |
| Conj, ColMajor, 1, Lower>::run), | |
| // array index: TR | (LO << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, RowMajor, false, Scalar, ColMajor, | |
| Conj, ColMajor, 1, Lower>::run), | |
| // array index: ADJ | (LO << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, RowMajor, Conj, Scalar, ColMajor, | |
| false, ColMajor, 1, Lower>::run), | |
| 0}; | |
| #endif | |
| const Scalar *a = reinterpret_cast<const Scalar *>(pa); | |
| Scalar *c = reinterpret_cast<Scalar *>(pc); | |
| Scalar alpha = *reinterpret_cast<const Scalar *>(palpha); | |
| Scalar beta = *reinterpret_cast<const Scalar *>(pbeta); | |
| EIGEN_BLAS_INT info = 0; | |
| if (UPLO(*uplo) == INVALID) | |
| info = 1; | |
| else if (OP(*op) == INVALID || (ISCOMPLEX && OP(*op) == ADJ)) | |
| info = 2; | |
| else if (*n < 0) | |
| info = 3; | |
| else if (*k < 0) | |
| info = 4; | |
| else if (*lda < std::max<EIGEN_BLAS_INT>(1, (OP(*op) == NOTR) ? *n : *k)) | |
| info = 7; | |
| else if (*ldc < std::max<EIGEN_BLAS_INT>(1, *n)) | |
| info = 10; | |
| if (info) return xerbla_(SCALAR_SUFFIX_UP "SYRK ", &info, kBlasNameLength); | |
| if (beta != Scalar(1)) { | |
| if (UPLO(*uplo) == UP) | |
| if (beta == Scalar(0)) | |
| matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero(); | |
| else | |
| matrix(c, *n, *n, *ldc).triangularView<Upper>() *= beta; | |
| else if (beta == Scalar(0)) | |
| matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero(); | |
| else | |
| matrix(c, *n, *n, *ldc).triangularView<Lower>() *= beta; | |
| } | |
| if (*n == 0 || *k == 0 || alpha == Scalar(0)) return; | |
| #if ISCOMPLEX | |
| // FIXME add support for symmetric complex matrix | |
| if (UPLO(*uplo) == UP) { | |
| if (OP(*op) == NOTR) | |
| matrix(c, *n, *n, *ldc).triangularView<Upper>() += | |
| alpha * matrix(a, *n, *k, *lda) * matrix(a, *n, *k, *lda).transpose(); | |
| else | |
| matrix(c, *n, *n, *ldc).triangularView<Upper>() += | |
| alpha * matrix(a, *k, *n, *lda).transpose() * matrix(a, *k, *n, *lda); | |
| } else { | |
| if (OP(*op) == NOTR) | |
| matrix(c, *n, *n, *ldc).triangularView<Lower>() += | |
| alpha * matrix(a, *n, *k, *lda) * matrix(a, *n, *k, *lda).transpose(); | |
| else | |
| matrix(c, *n, *n, *ldc).triangularView<Lower>() += | |
| alpha * matrix(a, *k, *n, *lda).transpose() * matrix(a, *k, *n, *lda); | |
| } | |
| #else | |
| Eigen::internal::gemm_blocking_space<ColMajor, Scalar, Scalar, Dynamic, Dynamic, Dynamic> 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); | |
| #endif | |
| } | |
| // c = alpha*a*b' + alpha*b*a' + beta*c for op = 'N'or'n' | |
| // c = alpha*a'*b + alpha*b'*a + beta*c for op = 'T'or't' | |
| EIGEN_BLAS_FUNC(syr2k) | |
| (const char *uplo, const char *op, const EIGEN_BLAS_INT *n, const EIGEN_BLAS_INT *k, const RealScalar *palpha, | |
| const RealScalar *pa, const EIGEN_BLAS_INT *lda, const RealScalar *pb, const EIGEN_BLAS_INT *ldb, | |
| const RealScalar *pbeta, RealScalar *pc, const EIGEN_BLAS_INT *ldc) { | |
| 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); | |
| EIGEN_BLAS_INT info = 0; | |
| if (UPLO(*uplo) == INVALID) | |
| info = 1; | |
| else if (OP(*op) == INVALID || (ISCOMPLEX && OP(*op) == ADJ)) | |
| info = 2; | |
| else if (*n < 0) | |
| info = 3; | |
| else if (*k < 0) | |
| info = 4; | |
| else if (*lda < std::max<EIGEN_BLAS_INT>(1, (OP(*op) == NOTR) ? *n : *k)) | |
| info = 7; | |
| else if (*ldb < std::max<EIGEN_BLAS_INT>(1, (OP(*op) == NOTR) ? *n : *k)) | |
| info = 9; | |
| else if (*ldc < std::max<EIGEN_BLAS_INT>(1, *n)) | |
| info = 12; | |
| if (info) return xerbla_(SCALAR_SUFFIX_UP "SYR2K", &info, kBlasNameLength); | |
| using Eigen::Lower; | |
| using Eigen::Upper; | |
| if (beta != Scalar(1)) { | |
| if (UPLO(*uplo) == UP) | |
| if (beta == Scalar(0)) | |
| matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero(); | |
| else | |
| matrix(c, *n, *n, *ldc).triangularView<Upper>() *= beta; | |
| else if (beta == Scalar(0)) | |
| matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero(); | |
| else | |
| matrix(c, *n, *n, *ldc).triangularView<Lower>() *= beta; | |
| } | |
| if (*k == 0) return; | |
| if (OP(*op) == NOTR) { | |
| if (UPLO(*uplo) == UP) { | |
| matrix(c, *n, *n, *ldc).triangularView<Upper>() += | |
| alpha * matrix(a, *n, *k, *lda) * matrix(b, *n, *k, *ldb).transpose() + | |
| alpha * matrix(b, *n, *k, *ldb) * matrix(a, *n, *k, *lda).transpose(); | |
| } else if (UPLO(*uplo) == LO) | |
| matrix(c, *n, *n, *ldc).triangularView<Lower>() += | |
| alpha * matrix(a, *n, *k, *lda) * matrix(b, *n, *k, *ldb).transpose() + | |
| alpha * matrix(b, *n, *k, *ldb) * matrix(a, *n, *k, *lda).transpose(); | |
| } else if (OP(*op) == TR || OP(*op) == ADJ) { | |
| if (UPLO(*uplo) == UP) | |
| matrix(c, *n, *n, *ldc).triangularView<Upper>() += | |
| alpha * matrix(a, *k, *n, *lda).transpose() * matrix(b, *k, *n, *ldb) + | |
| alpha * matrix(b, *k, *n, *ldb).transpose() * matrix(a, *k, *n, *lda); | |
| else if (UPLO(*uplo) == LO) | |
| matrix(c, *n, *n, *ldc).triangularView<Lower>() += | |
| alpha * matrix(a, *k, *n, *lda).transpose() * matrix(b, *k, *n, *ldb) + | |
| alpha * matrix(b, *k, *n, *ldb).transpose() * matrix(a, *k, *n, *lda); | |
| } | |
| } | |
| #if ISCOMPLEX | |
| // c = alpha*a*b + beta*c for side = 'L'or'l' | |
| // c = alpha*b*a + beta*c for side = 'R'or'r | |
| EIGEN_BLAS_FUNC(hemm) | |
| (const char *side, const char *uplo, 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 *ldb, | |
| const RealScalar *pbeta, RealScalar *pc, const EIGEN_BLAS_INT *ldc) { | |
| 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); | |
| EIGEN_BLAS_INT info = 0; | |
| if (SIDE(*side) == INVALID) | |
| info = 1; | |
| else if (UPLO(*uplo) == INVALID) | |
| info = 2; | |
| else if (*m < 0) | |
| info = 3; | |
| else if (*n < 0) | |
| info = 4; | |
| else if (*lda < std::max<EIGEN_BLAS_INT>(1, (SIDE(*side) == LEFT) ? *m : *n)) | |
| info = 7; | |
| else if (*ldb < std::max<EIGEN_BLAS_INT>(1, *m)) | |
| info = 9; | |
| else if (*ldc < std::max<EIGEN_BLAS_INT>(1, *m)) | |
| info = 12; | |
| if (info) return xerbla_(SCALAR_SUFFIX_UP "HEMM ", &info, kBlasNameLength); | |
| if (beta == Scalar(0)) | |
| matrix(c, *m, *n, *ldc).setZero(); | |
| else if (beta != Scalar(1)) | |
| matrix(c, *m, *n, *ldc) *= beta; | |
| if (*m == 0 || *n == 0 || alpha == Scalar(0)) return; | |
| using Eigen::ColMajor; | |
| using Eigen::DenseIndex; | |
| using Eigen::Dynamic; | |
| using Eigen::RowMajor; | |
| using Eigen::Upper; | |
| EIGEN_BLAS_INT size = (SIDE(*side) == LEFT) ? (*m) : (*n); | |
| Eigen::internal::gemm_blocking_space<ColMajor, Scalar, Scalar, Dynamic, Dynamic, Dynamic> blocking(*m, *n, size, 1, | |
| false); | |
| if (SIDE(*side) == LEFT) { | |
| if (UPLO(*uplo) == UP) | |
| Eigen::internal::product_selfadjoint_matrix<Scalar, DenseIndex, RowMajor, true, Conj, ColMajor, false, false, | |
| ColMajor, 1>::run(*m, *n, a, *lda, b, *ldb, c, 1, *ldc, alpha, | |
| blocking); | |
| else if (UPLO(*uplo) == LO) | |
| Eigen::internal::product_selfadjoint_matrix<Scalar, DenseIndex, ColMajor, true, false, ColMajor, false, false, | |
| ColMajor, 1>::run(*m, *n, a, *lda, b, *ldb, c, 1, *ldc, alpha, | |
| blocking); | |
| else | |
| return; | |
| } else if (SIDE(*side) == RIGHT) { | |
| if (UPLO(*uplo) == UP) | |
| matrix(c, *m, *n, *ldc) += | |
| alpha * matrix(b, *m, *n, *ldb) * | |
| matrix(a, *n, *n, *lda) | |
| .selfadjointView<Upper>(); /*internal::product_selfadjoint_matrix<Scalar,DenseIndex,ColMajor,false,false, | |
| RowMajor,true,Conj, ColMajor, 1> | |
| ::run(*m, *n, b, *ldb, a, *lda, c, 1, *ldc, alpha, blocking);*/ | |
| else if (UPLO(*uplo) == LO) | |
| Eigen::internal::product_selfadjoint_matrix<Scalar, DenseIndex, ColMajor, false, false, ColMajor, true, false, | |
| ColMajor, 1>::run(*m, *n, b, *ldb, a, *lda, c, 1, *ldc, alpha, | |
| blocking); | |
| else | |
| return; | |
| } else { | |
| return; | |
| } | |
| } | |
| // c = alpha*a*conj(a') + beta*c for op = 'N'or'n' | |
| // c = alpha*conj(a')*a + beta*c for op = 'C'or'c' | |
| EIGEN_BLAS_FUNC(herk) | |
| (const char *uplo, const char *op, const EIGEN_BLAS_INT *n, const EIGEN_BLAS_INT *k, const RealScalar *palpha, | |
| const RealScalar *pa, const EIGEN_BLAS_INT *lda, const RealScalar *pbeta, RealScalar *pc, const EIGEN_BLAS_INT *ldc) { | |
| 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<Scalar, Scalar> &); | |
| static const functype func[8] = { | |
| // array index: NOTR | (UP << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, ColMajor, false, Scalar, RowMajor, | |
| Conj, ColMajor, 1, Upper>::run), | |
| 0, | |
| // array index: ADJ | (UP << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, RowMajor, Conj, Scalar, ColMajor, | |
| false, ColMajor, 1, Upper>::run), | |
| 0, | |
| // array index: NOTR | (LO << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, ColMajor, false, Scalar, RowMajor, | |
| Conj, ColMajor, 1, Lower>::run), | |
| 0, | |
| // array index: ADJ | (LO << 2) | |
| (Eigen::internal::general_matrix_matrix_triangular_product<DenseIndex, Scalar, RowMajor, Conj, Scalar, ColMajor, | |
| false, ColMajor, 1, Lower>::run), | |
| 0}; | |
| const Scalar *a = reinterpret_cast<const Scalar *>(pa); | |
| Scalar *c = reinterpret_cast<Scalar *>(pc); | |
| RealScalar alpha = *palpha; | |
| RealScalar beta = *pbeta; | |
| EIGEN_BLAS_INT info = 0; | |
| if (UPLO(*uplo) == INVALID) | |
| info = 1; | |
| else if ((OP(*op) == INVALID) || (OP(*op) == TR)) | |
| info = 2; | |
| else if (*n < 0) | |
| info = 3; | |
| else if (*k < 0) | |
| info = 4; | |
| else if (*lda < std::max<EIGEN_BLAS_INT>(1, (OP(*op) == NOTR) ? *n : *k)) | |
| info = 7; | |
| else if (*ldc < std::max<EIGEN_BLAS_INT>(1, *n)) | |
| info = 10; | |
| if (info) return xerbla_(SCALAR_SUFFIX_UP "HERK ", &info, kBlasNameLength); | |
| int code = OP(*op) | (UPLO(*uplo) << 2); | |
| if (beta != RealScalar(1)) { | |
| if (UPLO(*uplo) == UP) | |
| if (beta == Scalar(0)) | |
| matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero(); | |
| else | |
| matrix(c, *n, *n, *ldc).triangularView<StrictlyUpper>() *= beta; | |
| else if (beta == Scalar(0)) | |
| matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero(); | |
| else | |
| matrix(c, *n, *n, *ldc).triangularView<StrictlyLower>() *= beta; | |
| if (beta != Scalar(0)) { | |
| matrix(c, *n, *n, *ldc).diagonal().real() *= beta; | |
| matrix(c, *n, *n, *ldc).diagonal().imag().setZero(); | |
| } | |
| } | |
| if (*k > 0 && alpha != RealScalar(0)) { | |
| Eigen::internal::gemm_blocking_space<ColMajor, Scalar, Scalar, Dynamic, Dynamic, Dynamic> 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(); | |
| } | |
| } | |
| // c = alpha*a*conj(b') + conj(alpha)*b*conj(a') + beta*c, for op = 'N'or'n' | |
| // c = alpha*conj(a')*b + conj(alpha)*conj(b')*a + beta*c, for op = 'C'or'c' | |
| EIGEN_BLAS_FUNC(her2k) | |
| (const char *uplo, const char *op, const EIGEN_BLAS_INT *n, const EIGEN_BLAS_INT *k, const RealScalar *palpha, | |
| const RealScalar *pa, const EIGEN_BLAS_INT *lda, const RealScalar *pb, const EIGEN_BLAS_INT *ldb, | |
| const RealScalar *pbeta, RealScalar *pc, const EIGEN_BLAS_INT *ldc) { | |
| 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); | |
| RealScalar beta = *pbeta; | |
| EIGEN_BLAS_INT info = 0; | |
| if (UPLO(*uplo) == INVALID) | |
| info = 1; | |
| else if ((OP(*op) == INVALID) || (OP(*op) == TR)) | |
| info = 2; | |
| else if (*n < 0) | |
| info = 3; | |
| else if (*k < 0) | |
| info = 4; | |
| else if (*lda < std::max<EIGEN_BLAS_INT>(1, (OP(*op) == NOTR) ? *n : *k)) | |
| info = 7; | |
| else if (*ldb < std::max<EIGEN_BLAS_INT>(1, (OP(*op) == NOTR) ? *n : *k)) | |
| info = 9; | |
| else if (*ldc < std::max<EIGEN_BLAS_INT>(1, *n)) | |
| info = 12; | |
| if (info) return xerbla_(SCALAR_SUFFIX_UP "HER2K", &info, kBlasNameLength); | |
| using Eigen::Lower; | |
| using Eigen::StrictlyLower; | |
| using Eigen::StrictlyUpper; | |
| using Eigen::Upper; | |
| if (beta != RealScalar(1)) { | |
| if (UPLO(*uplo) == UP) | |
| if (beta == Scalar(0)) | |
| matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero(); | |
| else | |
| matrix(c, *n, *n, *ldc).triangularView<StrictlyUpper>() *= beta; | |
| else if (beta == Scalar(0)) | |
| matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero(); | |
| else | |
| matrix(c, *n, *n, *ldc).triangularView<StrictlyLower>() *= beta; | |
| if (beta != Scalar(0)) { | |
| matrix(c, *n, *n, *ldc).diagonal().real() *= beta; | |
| matrix(c, *n, *n, *ldc).diagonal().imag().setZero(); | |
| } | |
| } else if (*k > 0 && alpha != Scalar(0)) | |
| matrix(c, *n, *n, *ldc).diagonal().imag().setZero(); | |
| if (*k == 0) return; | |
| if (OP(*op) == NOTR) { | |
| if (UPLO(*uplo) == UP) { | |
| matrix(c, *n, *n, *ldc).triangularView<Upper>() += | |
| alpha * matrix(a, *n, *k, *lda) * matrix(b, *n, *k, *ldb).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<Lower>() += | |
| alpha * matrix(a, *n, *k, *lda) * matrix(b, *n, *k, *ldb).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<Upper>() += | |
| alpha * matrix(a, *k, *n, *lda).adjoint() * matrix(b, *k, *n, *ldb) + | |
| 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<Lower>() += | |
| alpha * matrix(a, *k, *n, *lda).adjoint() * matrix(b, *k, *n, *ldb) + | |
| Eigen::numext::conj(alpha) * matrix(b, *k, *n, *ldb).adjoint() * matrix(a, *k, *n, *lda); | |
| } | |
| // The two products round differently, so force the documented real diagonal on exit. | |
| if (alpha != Scalar(0)) matrix(c, *n, *n, *ldc).diagonal().imag().setZero(); | |
| } | |
| #endif |