| /* |
| Copyright (c) 2011, Intel Corporation. All rights reserved. |
| |
| Redistribution and use in source and binary forms, with or without modification, |
| are permitted provided that the following conditions are met: |
| |
| * Redistributions of source code must retain the above copyright notice, this |
| list of conditions and the following disclaimer. |
| * Redistributions in binary form must reproduce the above copyright notice, |
| this list of conditions and the following disclaimer in the documentation |
| and/or other materials provided with the distribution. |
| * Neither the name of Intel Corporation nor the names of its contributors may |
| be used to endorse or promote products derived from this software without |
| specific prior written permission. |
| |
| THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND |
| ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED |
| WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE |
| DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR |
| ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES |
| (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; |
| LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON |
| ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT |
| (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS |
| SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
| |
| ******************************************************************************** |
| * Content : Eigen bindings to BLAS F77 |
| * Triangular matrix * matrix product functionality based on ?TRMM. |
| ******************************************************************************** |
| */ |
| |
| #ifndef EIGEN_TRIANGULAR_SOLVER_MATRIX_BLAS_H |
| #define EIGEN_TRIANGULAR_SOLVER_MATRIX_BLAS_H |
| |
| // IWYU pragma: private |
| #include "../InternalHeaderCheck.h" |
| |
| namespace Eigen { |
| |
| namespace internal { |
| |
| // implements LeftSide op(triangular)^-1 * general |
| #define EIGEN_BLAS_TRSM_L(EIGTYPE, BLASTYPE, BLASFUNC) \ |
| template <typename Index, int Mode, bool Conjugate, int TriStorageOrder> \ |
| struct triangular_solve_matrix<EIGTYPE, Index, OnTheLeft, Mode, Conjugate, TriStorageOrder, ColMajor, 1> { \ |
| enum { \ |
| IsLower = (Mode & Lower) == Lower, \ |
| IsUnitDiag = (Mode & UnitDiag) ? 1 : 0, \ |
| IsZeroDiag = (Mode & ZeroDiag) ? 1 : 0, \ |
| conjA = ((TriStorageOrder == ColMajor) && Conjugate) ? 1 : 0 \ |
| }; \ |
| static void run(Index size, Index otherSize, const EIGTYPE* _tri, Index triStride, EIGTYPE* _other, \ |
| Index otherIncr, Index otherStride, level3_blocking<EIGTYPE, EIGTYPE>& /*blocking*/) { \ |
| if (size == 0 || otherSize == 0) return; \ |
| EIGEN_ONLY_USED_FOR_DEBUG(otherIncr); \ |
| eigen_assert(otherIncr == 1); \ |
| BlasIndex m = convert_index<BlasIndex>(size), n = convert_index<BlasIndex>(otherSize), lda, ldb; \ |
| char side = 'L', uplo, diag = 'N', transa; \ |
| /* Set alpha_ */ \ |
| EIGTYPE alpha(1); \ |
| ldb = convert_index<BlasIndex>(otherStride); \ |
| \ |
| const EIGTYPE* a; \ |
| /* Set trans */ \ |
| transa = (TriStorageOrder == RowMajor) ? ((Conjugate) ? 'C' : 'T') : 'N'; \ |
| /* Set uplo */ \ |
| uplo = IsLower ? 'L' : 'U'; \ |
| if (TriStorageOrder == RowMajor) uplo = (uplo == 'L') ? 'U' : 'L'; \ |
| /* Set a, lda */ \ |
| typedef Matrix<EIGTYPE, Dynamic, Dynamic, TriStorageOrder> MatrixTri; \ |
| Map<const MatrixTri, 0, OuterStride<> > tri(_tri, size, size, OuterStride<>(triStride)); \ |
| MatrixTri a_tmp; \ |
| \ |
| if (conjA) { \ |
| a_tmp = tri.conjugate(); \ |
| a = a_tmp.data(); \ |
| lda = convert_index<BlasIndex>(a_tmp.outerStride()); \ |
| } else { \ |
| a = _tri; \ |
| lda = convert_index<BlasIndex>(triStride); \ |
| } \ |
| if (IsUnitDiag) diag = 'U'; \ |
| /* call ?trsm*/ \ |
| BLASFUNC(&side, &uplo, &transa, &diag, &m, &n, (const BLASTYPE*)&numext::real_ref(alpha), (const BLASTYPE*)a, \ |
| &lda, (BLASTYPE*)_other, &ldb); \ |
| } \ |
| }; |
| |
| #ifdef EIGEN_USE_MKL |
| EIGEN_BLAS_TRSM_L(double, double, dtrsm) |
| EIGEN_BLAS_TRSM_L(dcomplex, MKL_Complex16, ztrsm) |
| EIGEN_BLAS_TRSM_L(float, float, strsm) |
| EIGEN_BLAS_TRSM_L(scomplex, MKL_Complex8, ctrsm) |
| #else |
| EIGEN_BLAS_TRSM_L(double, double, dtrsm_) |
| EIGEN_BLAS_TRSM_L(dcomplex, double, ztrsm_) |
| EIGEN_BLAS_TRSM_L(float, float, strsm_) |
| EIGEN_BLAS_TRSM_L(scomplex, float, ctrsm_) |
| #endif |
| |
| // implements RightSide general * op(triangular)^-1 |
| #define EIGEN_BLAS_TRSM_R(EIGTYPE, BLASTYPE, BLASFUNC) \ |
| template <typename Index, int Mode, bool Conjugate, int TriStorageOrder> \ |
| struct triangular_solve_matrix<EIGTYPE, Index, OnTheRight, Mode, Conjugate, TriStorageOrder, ColMajor, 1> { \ |
| enum { \ |
| IsLower = (Mode & Lower) == Lower, \ |
| IsUnitDiag = (Mode & UnitDiag) ? 1 : 0, \ |
| IsZeroDiag = (Mode & ZeroDiag) ? 1 : 0, \ |
| conjA = ((TriStorageOrder == ColMajor) && Conjugate) ? 1 : 0 \ |
| }; \ |
| static void run(Index size, Index otherSize, const EIGTYPE* _tri, Index triStride, EIGTYPE* _other, \ |
| Index otherIncr, Index otherStride, level3_blocking<EIGTYPE, EIGTYPE>& /*blocking*/) { \ |
| if (size == 0 || otherSize == 0) return; \ |
| EIGEN_ONLY_USED_FOR_DEBUG(otherIncr); \ |
| eigen_assert(otherIncr == 1); \ |
| BlasIndex m = convert_index<BlasIndex>(otherSize), n = convert_index<BlasIndex>(size), lda, ldb; \ |
| char side = 'R', uplo, diag = 'N', transa; \ |
| /* Set alpha_ */ \ |
| EIGTYPE alpha(1); \ |
| ldb = convert_index<BlasIndex>(otherStride); \ |
| \ |
| const EIGTYPE* a; \ |
| /* Set trans */ \ |
| transa = (TriStorageOrder == RowMajor) ? ((Conjugate) ? 'C' : 'T') : 'N'; \ |
| /* Set uplo */ \ |
| uplo = IsLower ? 'L' : 'U'; \ |
| if (TriStorageOrder == RowMajor) uplo = (uplo == 'L') ? 'U' : 'L'; \ |
| /* Set a, lda */ \ |
| typedef Matrix<EIGTYPE, Dynamic, Dynamic, TriStorageOrder> MatrixTri; \ |
| Map<const MatrixTri, 0, OuterStride<> > tri(_tri, size, size, OuterStride<>(triStride)); \ |
| MatrixTri a_tmp; \ |
| \ |
| if (conjA) { \ |
| a_tmp = tri.conjugate(); \ |
| a = a_tmp.data(); \ |
| lda = convert_index<BlasIndex>(a_tmp.outerStride()); \ |
| } else { \ |
| a = _tri; \ |
| lda = convert_index<BlasIndex>(triStride); \ |
| } \ |
| if (IsUnitDiag) diag = 'U'; \ |
| /* call ?trsm*/ \ |
| BLASFUNC(&side, &uplo, &transa, &diag, &m, &n, (const BLASTYPE*)&numext::real_ref(alpha), (const BLASTYPE*)a, \ |
| &lda, (BLASTYPE*)_other, &ldb); \ |
| /*std::cout << "TRMS_L specialization!\n";*/ \ |
| } \ |
| }; |
| |
| #ifdef EIGEN_USE_MKL |
| EIGEN_BLAS_TRSM_R(double, double, dtrsm) |
| EIGEN_BLAS_TRSM_R(dcomplex, MKL_Complex16, ztrsm) |
| EIGEN_BLAS_TRSM_R(float, float, strsm) |
| EIGEN_BLAS_TRSM_R(scomplex, MKL_Complex8, ctrsm) |
| #else |
| EIGEN_BLAS_TRSM_R(double, double, dtrsm_) |
| EIGEN_BLAS_TRSM_R(dcomplex, double, ztrsm_) |
| EIGEN_BLAS_TRSM_R(float, float, strsm_) |
| EIGEN_BLAS_TRSM_R(scomplex, float, ctrsm_) |
| #endif |
| |
| } // end namespace internal |
| |
| } // end namespace Eigen |
| |
| #endif // EIGEN_TRIANGULAR_SOLVER_MATRIX_BLAS_H |
