From fec4c334bac76bfabd14168bf0ac668402f551a7 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 11 Apr 2016 16:04:09 +0200 Subject: Remove all references to MKL in BLAS wrappers. --- Eigen/Core | 16 +- .../products/GeneralMatrixMatrixTriangular_BLAS.h | 148 ++++++++++ .../products/GeneralMatrixMatrixTriangular_MKL.h | 148 ---------- Eigen/src/Core/products/GeneralMatrixMatrix_BLAS.h | 115 ++++++++ Eigen/src/Core/products/GeneralMatrixMatrix_MKL.h | 115 -------- Eigen/src/Core/products/GeneralMatrixVector_BLAS.h | 129 +++++++++ Eigen/src/Core/products/GeneralMatrixVector_MKL.h | 129 --------- .../Core/products/SelfadjointMatrixMatrix_BLAS.h | 275 +++++++++++++++++++ .../Core/products/SelfadjointMatrixMatrix_MKL.h | 275 ------------------- .../Core/products/SelfadjointMatrixVector_BLAS.h | 111 ++++++++ .../Core/products/SelfadjointMatrixVector_MKL.h | 111 -------- .../Core/products/TriangularMatrixMatrix_BLAS.h | 302 +++++++++++++++++++++ .../src/Core/products/TriangularMatrixMatrix_MKL.h | 302 --------------------- .../Core/products/TriangularMatrixVector_BLAS.h | 241 ++++++++++++++++ .../src/Core/products/TriangularMatrixVector_MKL.h | 241 ---------------- .../Core/products/TriangularSolverMatrix_BLAS.h | 151 +++++++++++ .../src/Core/products/TriangularSolverMatrix_MKL.h | 151 ----------- 17 files changed, 1480 insertions(+), 1480 deletions(-) create mode 100644 Eigen/src/Core/products/GeneralMatrixMatrixTriangular_BLAS.h delete mode 100644 Eigen/src/Core/products/GeneralMatrixMatrixTriangular_MKL.h create mode 100644 Eigen/src/Core/products/GeneralMatrixMatrix_BLAS.h delete mode 100644 Eigen/src/Core/products/GeneralMatrixMatrix_MKL.h create mode 100644 Eigen/src/Core/products/GeneralMatrixVector_BLAS.h delete mode 100644 Eigen/src/Core/products/GeneralMatrixVector_MKL.h create mode 100644 Eigen/src/Core/products/SelfadjointMatrixMatrix_BLAS.h delete mode 100644 Eigen/src/Core/products/SelfadjointMatrixMatrix_MKL.h create mode 100644 Eigen/src/Core/products/SelfadjointMatrixVector_BLAS.h delete mode 100644 Eigen/src/Core/products/SelfadjointMatrixVector_MKL.h create mode 100644 Eigen/src/Core/products/TriangularMatrixMatrix_BLAS.h delete mode 100644 Eigen/src/Core/products/TriangularMatrixMatrix_MKL.h create mode 100644 Eigen/src/Core/products/TriangularMatrixVector_BLAS.h delete mode 100644 Eigen/src/Core/products/TriangularMatrixVector_MKL.h create mode 100644 Eigen/src/Core/products/TriangularSolverMatrix_BLAS.h delete mode 100644 Eigen/src/Core/products/TriangularSolverMatrix_MKL.h (limited to 'Eigen') diff --git a/Eigen/Core b/Eigen/Core index 1e62f3ec1..30a572479 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -450,14 +450,14 @@ using std::ptrdiff_t; #include "src/Core/ArrayWrapper.h" #ifdef EIGEN_USE_BLAS -#include "src/Core/products/GeneralMatrixMatrix_MKL.h" -#include "src/Core/products/GeneralMatrixVector_MKL.h" -#include "src/Core/products/GeneralMatrixMatrixTriangular_MKL.h" -#include "src/Core/products/SelfadjointMatrixMatrix_MKL.h" -#include "src/Core/products/SelfadjointMatrixVector_MKL.h" -#include "src/Core/products/TriangularMatrixMatrix_MKL.h" -#include "src/Core/products/TriangularMatrixVector_MKL.h" -#include "src/Core/products/TriangularSolverMatrix_MKL.h" +#include "src/Core/products/GeneralMatrixMatrix_BLAS.h" +#include "src/Core/products/GeneralMatrixVector_BLAS.h" +#include "src/Core/products/GeneralMatrixMatrixTriangular_BLAS.h" +#include "src/Core/products/SelfadjointMatrixMatrix_BLAS.h" +#include "src/Core/products/SelfadjointMatrixVector_BLAS.h" +#include "src/Core/products/TriangularMatrixMatrix_BLAS.h" +#include "src/Core/products/TriangularMatrixVector_BLAS.h" +#include "src/Core/products/TriangularSolverMatrix_BLAS.h" #endif // EIGEN_USE_BLAS #ifdef EIGEN_USE_MKL_VML diff --git a/Eigen/src/Core/products/GeneralMatrixMatrixTriangular_BLAS.h b/Eigen/src/Core/products/GeneralMatrixMatrixTriangular_BLAS.h new file mode 100644 index 000000000..943d25bd1 --- /dev/null +++ b/Eigen/src/Core/products/GeneralMatrixMatrixTriangular_BLAS.h @@ -0,0 +1,148 @@ +/* + 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 + * Level 3 BLAS SYRK/HERK implementation. + ******************************************************************************** +*/ + +#ifndef EIGEN_GENERAL_MATRIX_MATRIX_TRIANGULAR_BLAS_H +#define EIGEN_GENERAL_MATRIX_MATRIX_TRIANGULAR_BLAS_H + +namespace Eigen { + +namespace internal { + +template +struct general_matrix_matrix_rankupdate : + general_matrix_matrix_triangular_product< + Index,Scalar,AStorageOrder,ConjugateA,Scalar,AStorageOrder,ConjugateA,ResStorageOrder,UpLo,BuiltIn> {}; + + +// try to go to BLAS specialization +#define EIGEN_BLAS_RANKUPDATE_SPECIALIZE(Scalar) \ +template \ +struct general_matrix_matrix_triangular_product { \ + static EIGEN_STRONG_INLINE void run(Index size, Index depth,const Scalar* lhs, Index lhsStride, \ + const Scalar* rhs, Index rhsStride, Scalar* res, Index resStride, Scalar alpha, level3_blocking& blocking) \ + { \ + if (lhs==rhs) { \ + general_matrix_matrix_rankupdate \ + ::run(size,depth,lhs,lhsStride,rhs,rhsStride,res,resStride,alpha,blocking); \ + } else { \ + general_matrix_matrix_triangular_product \ + ::run(size,depth,lhs,lhsStride,rhs,rhsStride,res,resStride,alpha,blocking); \ + } \ + } \ +}; + +EIGEN_BLAS_RANKUPDATE_SPECIALIZE(double) +EIGEN_BLAS_RANKUPDATE_SPECIALIZE(float) +// TODO handle complex cases +// EIGEN_BLAS_RANKUPDATE_SPECIALIZE(dcomplex) +// EIGEN_BLAS_RANKUPDATE_SPECIALIZE(scomplex) + +// SYRK for float/double +#define EIGEN_BLAS_RANKUPDATE_R(EIGTYPE, BLASTYPE, BLASFUNC) \ +template \ +struct general_matrix_matrix_rankupdate { \ + enum { \ + IsLower = (UpLo&Lower) == Lower, \ + LowUp = IsLower ? Lower : Upper, \ + conjA = ((AStorageOrder==ColMajor) && ConjugateA) ? 1 : 0 \ + }; \ + static EIGEN_STRONG_INLINE void run(Index size, Index depth,const EIGTYPE* lhs, Index lhsStride, \ + const EIGTYPE* /*rhs*/, Index /*rhsStride*/, EIGTYPE* res, Index resStride, EIGTYPE alpha, level3_blocking& /*blocking*/) \ + { \ + /* typedef Matrix MatrixRhs;*/ \ +\ + BlasIndex lda=convert_index(lhsStride), ldc=convert_index(resStride), n=convert_index(size), k=convert_index(depth); \ + char uplo=(IsLower) ? 'L' : 'U', trans=(AStorageOrder==RowMajor) ? 'T':'N'; \ + BLASTYPE alpha_, beta_; \ +\ +/* Set alpha_ & beta_ */ \ + assign_scalar_eig2mkl(alpha_, alpha); \ + assign_scalar_eig2mkl(beta_, EIGTYPE(1)); \ + BLASFUNC(&uplo, &trans, &n, &k, &alpha_, lhs, &lda, &beta_, res, &ldc); \ + } \ +}; + +// HERK for complex data +#define EIGEN_BLAS_RANKUPDATE_C(EIGTYPE, BLASTYPE, RTYPE, BLASFUNC) \ +template \ +struct general_matrix_matrix_rankupdate { \ + enum { \ + IsLower = (UpLo&Lower) == Lower, \ + LowUp = IsLower ? Lower : Upper, \ + conjA = (((AStorageOrder==ColMajor) && ConjugateA) || ((AStorageOrder==RowMajor) && !ConjugateA)) ? 1 : 0 \ + }; \ + static EIGEN_STRONG_INLINE void run(Index size, Index depth,const EIGTYPE* lhs, Index lhsStride, \ + const EIGTYPE* /*rhs*/, Index /*rhsStride*/, EIGTYPE* res, Index resStride, EIGTYPE alpha, level3_blocking& /*blocking*/) \ + { \ + typedef Matrix MatrixType; \ +\ + BlasIndex lda=convert_index(lhsStride), ldc=convert_index(resStride), n=convert_index(size), k=convert_index(depth); \ + char uplo=(IsLower) ? 'L' : 'U', trans=(AStorageOrder==RowMajor) ? 'C':'N'; \ + RTYPE alpha_, beta_; \ + const EIGTYPE* a_ptr; \ +\ +/* Set alpha_ & beta_ */ \ +/* assign_scalar_eig2mkl(alpha_, alpha); */\ +/* assign_scalar_eig2mkl(beta_, EIGTYPE(1));*/ \ + alpha_ = alpha.real(); \ + beta_ = 1.0; \ +/* Copy with conjugation in some cases*/ \ + MatrixType a; \ + if (conjA) { \ + Map > mapA(lhs,n,k,OuterStride<>(lhsStride)); \ + a = mapA.conjugate(); \ + lda = a.outerStride(); \ + a_ptr = a.data(); \ + } else a_ptr=lhs; \ + BLASFUNC(&uplo, &trans, &n, &k, &alpha_, (BLASTYPE*)a_ptr, &lda, &beta_, (BLASTYPE*)res, &ldc); \ + } \ +}; + + +EIGEN_BLAS_RANKUPDATE_R(double, double, dsyrk_) +EIGEN_BLAS_RANKUPDATE_R(float, float, ssyrk_) + +// TODO hanlde complex cases +// EIGEN_BLAS_RANKUPDATE_C(dcomplex, double, double, zherk_) +// EIGEN_BLAS_RANKUPDATE_C(scomplex, float, float, cherk_) + + +} // end namespace internal + +} // end namespace Eigen + +#endif // EIGEN_GENERAL_MATRIX_MATRIX_TRIANGULAR_BLAS_H diff --git a/Eigen/src/Core/products/GeneralMatrixMatrixTriangular_MKL.h b/Eigen/src/Core/products/GeneralMatrixMatrixTriangular_MKL.h deleted file mode 100644 index 6c835372c..000000000 --- a/Eigen/src/Core/products/GeneralMatrixMatrixTriangular_MKL.h +++ /dev/null @@ -1,148 +0,0 @@ -/* - 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 Intel(R) MKL - * Level 3 BLAS SYRK/HERK implementation. - ******************************************************************************** -*/ - -#ifndef EIGEN_GENERAL_MATRIX_MATRIX_TRIANGULAR_MKL_H -#define EIGEN_GENERAL_MATRIX_MATRIX_TRIANGULAR_MKL_H - -namespace Eigen { - -namespace internal { - -template -struct general_matrix_matrix_rankupdate : - general_matrix_matrix_triangular_product< - Index,Scalar,AStorageOrder,ConjugateA,Scalar,AStorageOrder,ConjugateA,ResStorageOrder,UpLo,BuiltIn> {}; - - -// try to go to BLAS specialization -#define EIGEN_MKL_RANKUPDATE_SPECIALIZE(Scalar) \ -template \ -struct general_matrix_matrix_triangular_product { \ - static EIGEN_STRONG_INLINE void run(Index size, Index depth,const Scalar* lhs, Index lhsStride, \ - const Scalar* rhs, Index rhsStride, Scalar* res, Index resStride, Scalar alpha, level3_blocking& blocking) \ - { \ - if (lhs==rhs) { \ - general_matrix_matrix_rankupdate \ - ::run(size,depth,lhs,lhsStride,rhs,rhsStride,res,resStride,alpha,blocking); \ - } else { \ - general_matrix_matrix_triangular_product \ - ::run(size,depth,lhs,lhsStride,rhs,rhsStride,res,resStride,alpha,blocking); \ - } \ - } \ -}; - -EIGEN_MKL_RANKUPDATE_SPECIALIZE(double) -EIGEN_MKL_RANKUPDATE_SPECIALIZE(float) -// TODO handle complex cases -// EIGEN_MKL_RANKUPDATE_SPECIALIZE(dcomplex) -// EIGEN_MKL_RANKUPDATE_SPECIALIZE(scomplex) - -// SYRK for float/double -#define EIGEN_MKL_RANKUPDATE_R(EIGTYPE, BLASTYPE, MKLFUNC) \ -template \ -struct general_matrix_matrix_rankupdate { \ - enum { \ - IsLower = (UpLo&Lower) == Lower, \ - LowUp = IsLower ? Lower : Upper, \ - conjA = ((AStorageOrder==ColMajor) && ConjugateA) ? 1 : 0 \ - }; \ - static EIGEN_STRONG_INLINE void run(Index size, Index depth,const EIGTYPE* lhs, Index lhsStride, \ - const EIGTYPE* /*rhs*/, Index /*rhsStride*/, EIGTYPE* res, Index resStride, EIGTYPE alpha, level3_blocking& /*blocking*/) \ - { \ - /* typedef Matrix MatrixRhs;*/ \ -\ - BlasIndex lda=convert_index(lhsStride), ldc=convert_index(resStride), n=convert_index(size), k=convert_index(depth); \ - char uplo=(IsLower) ? 'L' : 'U', trans=(AStorageOrder==RowMajor) ? 'T':'N'; \ - BLASTYPE alpha_, beta_; \ -\ -/* Set alpha_ & beta_ */ \ - assign_scalar_eig2mkl(alpha_, alpha); \ - assign_scalar_eig2mkl(beta_, EIGTYPE(1)); \ - MKLFUNC(&uplo, &trans, &n, &k, &alpha_, lhs, &lda, &beta_, res, &ldc); \ - } \ -}; - -// HERK for complex data -#define EIGEN_MKL_RANKUPDATE_C(EIGTYPE, BLASTYPE, RTYPE, MKLFUNC) \ -template \ -struct general_matrix_matrix_rankupdate { \ - enum { \ - IsLower = (UpLo&Lower) == Lower, \ - LowUp = IsLower ? Lower : Upper, \ - conjA = (((AStorageOrder==ColMajor) && ConjugateA) || ((AStorageOrder==RowMajor) && !ConjugateA)) ? 1 : 0 \ - }; \ - static EIGEN_STRONG_INLINE void run(Index size, Index depth,const EIGTYPE* lhs, Index lhsStride, \ - const EIGTYPE* /*rhs*/, Index /*rhsStride*/, EIGTYPE* res, Index resStride, EIGTYPE alpha, level3_blocking& /*blocking*/) \ - { \ - typedef Matrix MatrixType; \ -\ - BlasIndex lda=convert_index(lhsStride), ldc=convert_index(resStride), n=convert_index(size), k=convert_index(depth); \ - char uplo=(IsLower) ? 'L' : 'U', trans=(AStorageOrder==RowMajor) ? 'C':'N'; \ - RTYPE alpha_, beta_; \ - const EIGTYPE* a_ptr; \ -\ -/* Set alpha_ & beta_ */ \ -/* assign_scalar_eig2mkl(alpha_, alpha); */\ -/* assign_scalar_eig2mkl(beta_, EIGTYPE(1));*/ \ - alpha_ = alpha.real(); \ - beta_ = 1.0; \ -/* Copy with conjugation in some cases*/ \ - MatrixType a; \ - if (conjA) { \ - Map > mapA(lhs,n,k,OuterStride<>(lhsStride)); \ - a = mapA.conjugate(); \ - lda = a.outerStride(); \ - a_ptr = a.data(); \ - } else a_ptr=lhs; \ - MKLFUNC(&uplo, &trans, &n, &k, &alpha_, (BLASTYPE*)a_ptr, &lda, &beta_, (BLASTYPE*)res, &ldc); \ - } \ -}; - - -EIGEN_MKL_RANKUPDATE_R(double, double, dsyrk_) -EIGEN_MKL_RANKUPDATE_R(float, float, ssyrk_) - -// TODO hanlde complex cases -// EIGEN_MKL_RANKUPDATE_C(dcomplex, double, double, zherk_) -// EIGEN_MKL_RANKUPDATE_C(scomplex, float, float, cherk_) - - -} // end namespace internal - -} // end namespace Eigen - -#endif // EIGEN_GENERAL_MATRIX_MATRIX_TRIANGULAR_MKL_H diff --git a/Eigen/src/Core/products/GeneralMatrixMatrix_BLAS.h b/Eigen/src/Core/products/GeneralMatrixMatrix_BLAS.h new file mode 100644 index 000000000..7a3bdbf20 --- /dev/null +++ b/Eigen/src/Core/products/GeneralMatrixMatrix_BLAS.h @@ -0,0 +1,115 @@ +/* + 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 + * General matrix-matrix product functionality based on ?GEMM. + ******************************************************************************** +*/ + +#ifndef EIGEN_GENERAL_MATRIX_MATRIX_BLAS_H +#define EIGEN_GENERAL_MATRIX_MATRIX_BLAS_H + +namespace Eigen { + +namespace internal { + +/********************************************************************** +* This file implements general matrix-matrix multiplication using BLAS +* gemm function via partial specialization of +* general_matrix_matrix_product::run(..) method for float, double, +* std::complex and std::complex types +**********************************************************************/ + +// gemm specialization + +#define GEMM_SPECIALIZATION(EIGTYPE, EIGPREFIX, BLASTYPE, BLASPREFIX) \ +template< \ + typename Index, \ + int LhsStorageOrder, bool ConjugateLhs, \ + int RhsStorageOrder, bool ConjugateRhs> \ +struct general_matrix_matrix_product \ +{ \ +typedef gebp_traits Traits; \ +\ +static void run(Index rows, Index cols, Index depth, \ + const EIGTYPE* _lhs, Index lhsStride, \ + const EIGTYPE* _rhs, Index rhsStride, \ + EIGTYPE* res, Index resStride, \ + EIGTYPE alpha, \ + level3_blocking& /*blocking*/, \ + GemmParallelInfo* /*info = 0*/) \ +{ \ + using std::conj; \ +\ + char transa, transb; \ + BlasIndex m, n, k, lda, ldb, ldc; \ + const EIGTYPE *a, *b; \ + EIGTYPE beta(1); \ + MatrixX##EIGPREFIX a_tmp, b_tmp; \ +\ +/* Set transpose options */ \ + transa = (LhsStorageOrder==RowMajor) ? ((ConjugateLhs) ? 'C' : 'T') : 'N'; \ + transb = (RhsStorageOrder==RowMajor) ? ((ConjugateRhs) ? 'C' : 'T') : 'N'; \ +\ +/* Set m, n, k */ \ + m = convert_index(rows); \ + n = convert_index(cols); \ + k = convert_index(depth); \ +\ +/* Set lda, ldb, ldc */ \ + lda = convert_index(lhsStride); \ + ldb = convert_index(rhsStride); \ + ldc = convert_index(resStride); \ +\ +/* Set a, b, c */ \ + if ((LhsStorageOrder==ColMajor) && (ConjugateLhs)) { \ + Map > lhs(_lhs,m,k,OuterStride<>(lhsStride)); \ + a_tmp = lhs.conjugate(); \ + a = a_tmp.data(); \ + lda = convert_index(a_tmp.outerStride()); \ + } else a = _lhs; \ +\ + if ((RhsStorageOrder==ColMajor) && (ConjugateRhs)) { \ + Map > rhs(_rhs,k,n,OuterStride<>(rhsStride)); \ + b_tmp = rhs.conjugate(); \ + b = b_tmp.data(); \ + ldb = convert_index(b_tmp.outerStride()); \ + } else b = _rhs; \ +\ + BLASPREFIX##gemm_(&transa, &transb, &m, &n, &k, &numext::real_ref(alpha), (const BLASTYPE*)a, &lda, (const BLASTYPE*)b, &ldb, &numext::real_ref(beta), (BLASTYPE*)res, &ldc); \ +}}; + +GEMM_SPECIALIZATION(double, d, double, d) +GEMM_SPECIALIZATION(float, f, float, s) +GEMM_SPECIALIZATION(dcomplex, cd, double, z) +GEMM_SPECIALIZATION(scomplex, cf, float, c) + +} // end namespase internal + +} // end namespace Eigen + +#endif // EIGEN_GENERAL_MATRIX_MATRIX_BLAS_H diff --git a/Eigen/src/Core/products/GeneralMatrixMatrix_MKL.h b/Eigen/src/Core/products/GeneralMatrixMatrix_MKL.h deleted file mode 100644 index 299faf2f2..000000000 --- a/Eigen/src/Core/products/GeneralMatrixMatrix_MKL.h +++ /dev/null @@ -1,115 +0,0 @@ -/* - 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 Intel(R) MKL - * General matrix-matrix product functionality based on ?GEMM. - ******************************************************************************** -*/ - -#ifndef EIGEN_GENERAL_MATRIX_MATRIX_MKL_H -#define EIGEN_GENERAL_MATRIX_MATRIX_MKL_H - -namespace Eigen { - -namespace internal { - -/********************************************************************** -* This file implements general matrix-matrix multiplication using BLAS -* gemm function via partial specialization of -* general_matrix_matrix_product::run(..) method for float, double, -* std::complex and std::complex types -**********************************************************************/ - -// gemm specialization - -#define GEMM_SPECIALIZATION(EIGTYPE, EIGPREFIX, BLASTYPE, MKLPREFIX) \ -template< \ - typename Index, \ - int LhsStorageOrder, bool ConjugateLhs, \ - int RhsStorageOrder, bool ConjugateRhs> \ -struct general_matrix_matrix_product \ -{ \ -typedef gebp_traits Traits; \ -\ -static void run(Index rows, Index cols, Index depth, \ - const EIGTYPE* _lhs, Index lhsStride, \ - const EIGTYPE* _rhs, Index rhsStride, \ - EIGTYPE* res, Index resStride, \ - EIGTYPE alpha, \ - level3_blocking& /*blocking*/, \ - GemmParallelInfo* /*info = 0*/) \ -{ \ - using std::conj; \ -\ - char transa, transb; \ - BlasIndex m, n, k, lda, ldb, ldc; \ - const EIGTYPE *a, *b; \ - EIGTYPE beta(1); \ - MatrixX##EIGPREFIX a_tmp, b_tmp; \ -\ -/* Set transpose options */ \ - transa = (LhsStorageOrder==RowMajor) ? ((ConjugateLhs) ? 'C' : 'T') : 'N'; \ - transb = (RhsStorageOrder==RowMajor) ? ((ConjugateRhs) ? 'C' : 'T') : 'N'; \ -\ -/* Set m, n, k */ \ - m = convert_index(rows); \ - n = convert_index(cols); \ - k = convert_index(depth); \ -\ -/* Set lda, ldb, ldc */ \ - lda = convert_index(lhsStride); \ - ldb = convert_index(rhsStride); \ - ldc = convert_index(resStride); \ -\ -/* Set a, b, c */ \ - if ((LhsStorageOrder==ColMajor) && (ConjugateLhs)) { \ - Map > lhs(_lhs,m,k,OuterStride<>(lhsStride)); \ - a_tmp = lhs.conjugate(); \ - a = a_tmp.data(); \ - lda = convert_index(a_tmp.outerStride()); \ - } else a = _lhs; \ -\ - if ((RhsStorageOrder==ColMajor) && (ConjugateRhs)) { \ - Map > rhs(_rhs,k,n,OuterStride<>(rhsStride)); \ - b_tmp = rhs.conjugate(); \ - b = b_tmp.data(); \ - ldb = convert_index(b_tmp.outerStride()); \ - } else b = _rhs; \ -\ - MKLPREFIX##gemm_(&transa, &transb, &m, &n, &k, &numext::real_ref(alpha), (const BLASTYPE*)a, &lda, (const BLASTYPE*)b, &ldb, &numext::real_ref(beta), (BLASTYPE*)res, &ldc); \ -}}; - -GEMM_SPECIALIZATION(double, d, double, d) -GEMM_SPECIALIZATION(float, f, float, s) -GEMM_SPECIALIZATION(dcomplex, cd, double, z) -GEMM_SPECIALIZATION(scomplex, cf, float, c) - -} // end namespase internal - -} // end namespace Eigen - -#endif // EIGEN_GENERAL_MATRIX_MATRIX_MKL_H diff --git a/Eigen/src/Core/products/GeneralMatrixVector_BLAS.h b/Eigen/src/Core/products/GeneralMatrixVector_BLAS.h new file mode 100644 index 000000000..e3a5d5892 --- /dev/null +++ b/Eigen/src/Core/products/GeneralMatrixVector_BLAS.h @@ -0,0 +1,129 @@ +/* + 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 + * General matrix-vector product functionality based on ?GEMV. + ******************************************************************************** +*/ + +#ifndef EIGEN_GENERAL_MATRIX_VECTOR_BLAS_H +#define EIGEN_GENERAL_MATRIX_VECTOR_BLAS_H + +namespace Eigen { + +namespace internal { + +/********************************************************************** +* This file implements general matrix-vector multiplication using BLAS +* gemv function via partial specialization of +* general_matrix_vector_product::run(..) method for float, double, +* std::complex and std::complex types +**********************************************************************/ + +// gemv specialization + +template +struct general_matrix_vector_product_gemv; + +#define EIGEN_BLAS_GEMV_SPECIALIZE(Scalar) \ +template \ +struct general_matrix_vector_product,ColMajor,ConjugateLhs,Scalar,const_blas_data_mapper,ConjugateRhs,Specialized> { \ +static void run( \ + Index rows, Index cols, \ + const const_blas_data_mapper &lhs, \ + const const_blas_data_mapper &rhs, \ + Scalar* res, Index resIncr, Scalar alpha) \ +{ \ + if (ConjugateLhs) { \ + general_matrix_vector_product,ColMajor,ConjugateLhs,Scalar,const_blas_data_mapper,ConjugateRhs,BuiltIn>::run( \ + rows, cols, lhs, rhs, res, resIncr, alpha); \ + } else { \ + general_matrix_vector_product_gemv::run( \ + rows, cols, lhs.data(), lhs.stride(), rhs.data(), rhs.stride(), res, resIncr, alpha); \ + } \ +} \ +}; \ +template \ +struct general_matrix_vector_product,RowMajor,ConjugateLhs,Scalar,const_blas_data_mapper,ConjugateRhs,Specialized> { \ +static void run( \ + Index rows, Index cols, \ + const const_blas_data_mapper &lhs, \ + const const_blas_data_mapper &rhs, \ + Scalar* res, Index resIncr, Scalar alpha) \ +{ \ + general_matrix_vector_product_gemv::run( \ + rows, cols, lhs.data(), lhs.stride(), rhs.data(), rhs.stride(), res, resIncr, alpha); \ +} \ +}; \ + +EIGEN_BLAS_GEMV_SPECIALIZE(double) +EIGEN_BLAS_GEMV_SPECIALIZE(float) +EIGEN_BLAS_GEMV_SPECIALIZE(dcomplex) +EIGEN_BLAS_GEMV_SPECIALIZE(scomplex) + +#define EIGEN_BLAS_GEMV_SPECIALIZATION(EIGTYPE,BLASTYPE,BLASPREFIX) \ +template \ +struct general_matrix_vector_product_gemv \ +{ \ +typedef Matrix GEMVVector;\ +\ +static void run( \ + Index rows, Index cols, \ + const EIGTYPE* lhs, Index lhsStride, \ + const EIGTYPE* rhs, Index rhsIncr, \ + EIGTYPE* res, Index resIncr, EIGTYPE alpha) \ +{ \ + BlasIndex m=convert_index(rows), n=convert_index(cols), \ + lda=convert_index(lhsStride), incx=convert_index(rhsIncr), incy=convert_index(resIncr); \ + const EIGTYPE beta(1); \ + const EIGTYPE *x_ptr; \ + char trans=(LhsStorageOrder==ColMajor) ? 'N' : (ConjugateLhs) ? 'C' : 'T'; \ + if (LhsStorageOrder==RowMajor) { \ + m = convert_index(cols); \ + n = convert_index(rows); \ + }\ + GEMVVector x_tmp; \ + if (ConjugateRhs) { \ + Map > map_x(rhs,cols,1,InnerStride<>(incx)); \ + x_tmp=map_x.conjugate(); \ + x_ptr=x_tmp.data(); \ + incx=1; \ + } else x_ptr=rhs; \ + BLASPREFIX##gemv_(&trans, &m, &n, &numext::real_ref(alpha), (const BLASTYPE*)lhs, &lda, (const BLASTYPE*)x_ptr, &incx, &numext::real_ref(beta), (BLASTYPE*)res, &incy); \ +}\ +}; + +EIGEN_BLAS_GEMV_SPECIALIZATION(double, double, d) +EIGEN_BLAS_GEMV_SPECIALIZATION(float, float, s) +EIGEN_BLAS_GEMV_SPECIALIZATION(dcomplex, double, z) +EIGEN_BLAS_GEMV_SPECIALIZATION(scomplex, float, c) + +} // end namespase internal + +} // end namespace Eigen + +#endif // EIGEN_GENERAL_MATRIX_VECTOR_BLAS_H diff --git a/Eigen/src/Core/products/GeneralMatrixVector_MKL.h b/Eigen/src/Core/products/GeneralMatrixVector_MKL.h deleted file mode 100644 index c447c4aed..000000000 --- a/Eigen/src/Core/products/GeneralMatrixVector_MKL.h +++ /dev/null @@ -1,129 +0,0 @@ -/* - 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 Intel(R) MKL - * General matrix-vector product functionality based on ?GEMV. - ******************************************************************************** -*/ - -#ifndef EIGEN_GENERAL_MATRIX_VECTOR_MKL_H -#define EIGEN_GENERAL_MATRIX_VECTOR_MKL_H - -namespace Eigen { - -namespace internal { - -/********************************************************************** -* This file implements general matrix-vector multiplication using BLAS -* gemv function via partial specialization of -* general_matrix_vector_product::run(..) method for float, double, -* std::complex and std::complex types -**********************************************************************/ - -// gemv specialization - -template -struct general_matrix_vector_product_gemv; - -#define EIGEN_MKL_GEMV_SPECIALIZE(Scalar) \ -template \ -struct general_matrix_vector_product,ColMajor,ConjugateLhs,Scalar,const_blas_data_mapper,ConjugateRhs,Specialized> { \ -static void run( \ - Index rows, Index cols, \ - const const_blas_data_mapper &lhs, \ - const const_blas_data_mapper &rhs, \ - Scalar* res, Index resIncr, Scalar alpha) \ -{ \ - if (ConjugateLhs) { \ - general_matrix_vector_product,ColMajor,ConjugateLhs,Scalar,const_blas_data_mapper,ConjugateRhs,BuiltIn>::run( \ - rows, cols, lhs, rhs, res, resIncr, alpha); \ - } else { \ - general_matrix_vector_product_gemv::run( \ - rows, cols, lhs.data(), lhs.stride(), rhs.data(), rhs.stride(), res, resIncr, alpha); \ - } \ -} \ -}; \ -template \ -struct general_matrix_vector_product,RowMajor,ConjugateLhs,Scalar,const_blas_data_mapper,ConjugateRhs,Specialized> { \ -static void run( \ - Index rows, Index cols, \ - const const_blas_data_mapper &lhs, \ - const const_blas_data_mapper &rhs, \ - Scalar* res, Index resIncr, Scalar alpha) \ -{ \ - general_matrix_vector_product_gemv::run( \ - rows, cols, lhs.data(), lhs.stride(), rhs.data(), rhs.stride(), res, resIncr, alpha); \ -} \ -}; \ - -EIGEN_MKL_GEMV_SPECIALIZE(double) -EIGEN_MKL_GEMV_SPECIALIZE(float) -EIGEN_MKL_GEMV_SPECIALIZE(dcomplex) -EIGEN_MKL_GEMV_SPECIALIZE(scomplex) - -#define EIGEN_MKL_GEMV_SPECIALIZATION(EIGTYPE,BLASTYPE,MKLPREFIX) \ -template \ -struct general_matrix_vector_product_gemv \ -{ \ -typedef Matrix GEMVVector;\ -\ -static void run( \ - Index rows, Index cols, \ - const EIGTYPE* lhs, Index lhsStride, \ - const EIGTYPE* rhs, Index rhsIncr, \ - EIGTYPE* res, Index resIncr, EIGTYPE alpha) \ -{ \ - BlasIndex m=convert_index(rows), n=convert_index(cols), \ - lda=convert_index(lhsStride), incx=convert_index(rhsIncr), incy=convert_index(resIncr); \ - const EIGTYPE beta(1); \ - const EIGTYPE *x_ptr; \ - char trans=(LhsStorageOrder==ColMajor) ? 'N' : (ConjugateLhs) ? 'C' : 'T'; \ - if (LhsStorageOrder==RowMajor) { \ - m = convert_index(cols); \ - n = convert_index(rows); \ - }\ - GEMVVector x_tmp; \ - if (ConjugateRhs) { \ - Map > map_x(rhs,cols,1,InnerStride<>(incx)); \ - x_tmp=map_x.conjugate(); \ - x_ptr=x_tmp.data(); \ - incx=1; \ - } else x_ptr=rhs; \ - MKLPREFIX##gemv_(&trans, &m, &n, &numext::real_ref(alpha), (const BLASTYPE*)lhs, &lda, (const BLASTYPE*)x_ptr, &incx, &numext::real_ref(beta), (BLASTYPE*)res, &incy); \ -}\ -}; - -EIGEN_MKL_GEMV_SPECIALIZATION(double, double, d) -EIGEN_MKL_GEMV_SPECIALIZATION(float, float, s) -EIGEN_MKL_GEMV_SPECIALIZATION(dcomplex, double, z) -EIGEN_MKL_GEMV_SPECIALIZATION(scomplex, float, c) - -} // end namespase internal - -} // end namespace Eigen - -#endif // EIGEN_GENERAL_MATRIX_VECTOR_MKL_H diff --git a/Eigen/src/Core/products/SelfadjointMatrixMatrix_BLAS.h b/Eigen/src/Core/products/SelfadjointMatrixMatrix_BLAS.h new file mode 100644 index 000000000..c3e37b1e0 --- /dev/null +++ b/Eigen/src/Core/products/SelfadjointMatrixMatrix_BLAS.h @@ -0,0 +1,275 @@ +/* + 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 + * Self adjoint matrix * matrix product functionality based on ?SYMM/?HEMM. + ******************************************************************************** +*/ + +#ifndef EIGEN_SELFADJOINT_MATRIX_MATRIX_BLAS_H +#define EIGEN_SELFADJOINT_MATRIX_MATRIX_BLAS_H + +namespace Eigen { + +namespace internal { + + +/* Optimized selfadjoint matrix * matrix (?SYMM/?HEMM) product */ + +#define EIGEN_BLAS_SYMM_L(EIGTYPE, BLASTYPE, EIGPREFIX, BLASPREFIX) \ +template \ +struct product_selfadjoint_matrix \ +{\ +\ + static void run( \ + Index rows, Index cols, \ + const EIGTYPE* _lhs, Index lhsStride, \ + const EIGTYPE* _rhs, Index rhsStride, \ + EIGTYPE* res, Index resStride, \ + EIGTYPE alpha, level3_blocking& /*blocking*/) \ + { \ + char side='L', uplo='L'; \ + BlasIndex m, n, lda, ldb, ldc; \ + const EIGTYPE *a, *b; \ + EIGTYPE beta(1); \ + MatrixX##EIGPREFIX b_tmp; \ +\ +/* Set transpose options */ \ +/* Set m, n, k */ \ + m = convert_index(rows); \ + n = convert_index(cols); \ +\ +/* Set lda, ldb, ldc */ \ + lda = convert_index(lhsStride); \ + ldb = convert_index(rhsStride); \ + ldc = convert_index(resStride); \ +\ +/* Set a, b, c */ \ + if (LhsStorageOrder==RowMajor) uplo='U'; \ + a = _lhs; \ +\ + if (RhsStorageOrder==RowMajor) { \ + Map > rhs(_rhs,n,m,OuterStride<>(rhsStride)); \ + b_tmp = rhs.adjoint(); \ + b = b_tmp.data(); \ + ldb = convert_index(b_tmp.outerStride()); \ + } else b = _rhs; \ +\ + BLASPREFIX##symm_(&side, &uplo, &m, &n, &numext::real_ref(alpha), (const BLASTYPE*)a, &lda, (const BLASTYPE*)b, &ldb, &numext::real_ref(beta), (BLASTYPE*)res, &ldc); \ +\ + } \ +}; + + +#define EIGEN_BLAS_HEMM_L(EIGTYPE, BLASTYPE, EIGPREFIX, BLASPREFIX) \ +template \ +struct product_selfadjoint_matrix \ +{\ + static void run( \ + Index rows, Index cols, \ + const EIGTYPE* _lhs, Index lhsStride, \ + const EIGTYPE* _rhs, Index rhsStride, \ + EIGTYPE* res, Index resStride, \ + EIGTYPE alpha, level3_blocking& /*blocking*/) \ + { \ + char side='L', uplo='L'; \ + BlasIndex m, n, lda, ldb, ldc; \ + const EIGTYPE *a, *b; \ + EIGTYPE beta(1); \ + MatrixX##EIGPREFIX b_tmp; \ + Matrix a_tmp; \ +\ +/* Set transpose options */ \ +/* Set m, n, k */ \ + m = convert_index(rows); \ + n = convert_index(cols); \ +\ +/* Set lda, ldb, ldc */ \ + lda = convert_index(lhsStride); \ + ldb = convert_index(rhsStride); \ + ldc = convert_index(resStride); \ +\ +/* Set a, b, c */ \ + if (((LhsStorageOrder==ColMajor) && ConjugateLhs) || ((LhsStorageOrder==RowMajor) && (!ConjugateLhs))) { \ + Map, 0, OuterStride<> > lhs(_lhs,m,m,OuterStride<>(lhsStride)); \ + a_tmp = lhs.conjugate(); \ + a = a_tmp.data(); \ + lda = a_tmp.outerStride(); \ + } else a = _lhs; \ + if (LhsStorageOrder==RowMajor) uplo='U'; \ +\ + if (RhsStorageOrder==ColMajor && (!ConjugateRhs)) { \ + b = _rhs; } \ + else { \ + if (RhsStorageOrder==ColMajor && ConjugateRhs) { \ + Map > rhs(_rhs,m,n,OuterStride<>(rhsStride)); \ + b_tmp = rhs.conjugate(); \ + } else \ + if (ConjugateRhs) { \ + Map > rhs(_rhs,n,m,OuterStride<>(rhsStride)); \ + b_tmp = rhs.adjoint(); \ + } else { \ + Map > rhs(_rhs,n,m,OuterStride<>(rhsStride)); \ + b_tmp = rhs.transpose(); \ + } \ + b = b_tmp.data(); \ + ldb = convert_index(b_tmp.outerStride()); \ + } \ +\ + BLASPREFIX##hemm_(&side, &uplo, &m, &n, &numext::real_ref(alpha), (const BLASTYPE*)a, &lda, (const BLASTYPE*)b, &ldb, &numext::real_ref(beta), (BLASTYPE*)res, &ldc); \ +\ + } \ +}; + +EIGEN_BLAS_SYMM_L(double, double, d, d) +EIGEN_BLAS_SYMM_L(float, float, f, s) +EIGEN_BLAS_HEMM_L(dcomplex, double, cd, z) +EIGEN_BLAS_HEMM_L(scomplex, float, cf, c) + + +/* Optimized matrix * selfadjoint matrix (?SYMM/?HEMM) product */ + +#define EIGEN_BLAS_SYMM_R(EIGTYPE, BLASTYPE, EIGPREFIX, BLASPREFIX) \ +template \ +struct product_selfadjoint_matrix \ +{\ +\ + static void run( \ + Index rows, Index cols, \ + const EIGTYPE* _lhs, Index lhsStride, \ + const EIGTYPE* _rhs, Index rhsStride, \ + EIGTYPE* res, Index resStride, \ + EIGTYPE alpha, level3_blocking& /*blocking*/) \ + { \ + char side='R', uplo='L'; \ + BlasIndex m, n, lda, ldb, ldc; \ + const EIGTYPE *a, *b; \ + EIGTYPE beta(1); \ + MatrixX##EIGPREFIX b_tmp; \ +\ +/* Set m, n, k */ \ + m = convert_index(rows); \ + n = convert_index(cols); \ +\ +/* Set lda, ldb, ldc */ \ + lda = convert_index(rhsStride); \ + ldb = convert_index(lhsStride); \ + ldc = convert_index(resStride); \ +\ +/* Set a, b, c */ \ + if (RhsStorageOrder==RowMajor) uplo='U'; \ + a = _rhs; \ +\ + if (LhsStorageOrder==RowMajor) { \ + Map > lhs(_lhs,n,m,OuterStride<>(rhsStride)); \ + b_tmp = lhs.adjoint(); \ + b = b_tmp.data(); \ + ldb = convert_index(b_tmp.outerStride()); \ + } else b = _lhs; \ +\ + BLASPREFIX##symm_(&side, &uplo, &m, &n, &numext::real_ref(alpha), (const BLASTYPE*)a, &lda, (const BLASTYPE*)b, &ldb, &numext::real_ref(beta), (BLASTYPE*)res, &ldc); \ +\ + } \ +}; + + +#define EIGEN_BLAS_HEMM_R(EIGTYPE, BLASTYPE, EIGPREFIX, BLASPREFIX) \ +template \ +struct product_selfadjoint_matrix \ +{\ + static void run( \ + Index rows, Index cols, \ + const EIGTYPE* _lhs, Index lhsStride, \ + const EIGTYPE* _rhs, Index rhsStride, \ + EIGTYPE* res, Index resStride, \ + EIGTYPE alpha, level3_blocking& /*blocking*/) \ + { \ + char side='R', uplo='L'; \ + BlasIndex m, n, lda, ldb, ldc; \ + const EIGTYPE *a, *b; \ + EIGTYPE beta(1); \ + MatrixX##EIGPREFIX b_tmp; \ + Matrix a_tmp; \ +\ +/* Set m, n, k */ \ + m = convert_index(rows); \ + n = convert_index(cols); \ +\ +/* Set lda, ldb, ldc */ \ + lda = convert_index(rhsStride); \ + ldb = convert_index(lhsStride); \ + ldc = convert_index(resStride); \ +\ +/* Set a, b, c */ \ + if (((RhsStorageOrder==ColMajor) && ConjugateRhs) || ((RhsStorageOrder==RowMajor) && (!ConjugateRhs))) { \ + Map, 0, OuterStride<> > rhs(_rhs,n,n,OuterStride<>(rhsStride)); \ + a_tmp = rhs.conjugate(); \ + a = a_tmp.data(); \ + lda = convert_index(a_tmp.outerStride()); \ + } else a = _rhs; \ + if (RhsStorageOrder==RowMajor) uplo='U'; \ +\ + if (LhsStorageOrder==ColMajor && (!ConjugateLhs)) { \ + b = _lhs; } \ + else { \ + if (LhsStorageOrder==ColMajor && ConjugateLhs) { \ + Map > lhs(_lhs,m,n,OuterStride<>(lhsStride)); \ + b_tmp = lhs.conjugate(); \ + } else \ + if (ConjugateLhs) { \ + Map > lhs(_lhs,n,m,OuterStride<>(lhsStride)); \ + b_tmp = lhs.adjoint(); \ + } else { \ + Map > lhs(_lhs,n,m,OuterStride<>(lhsStride)); \ + b_tmp = lhs.transpose(); \ + } \ + b = b_tmp.data(); \ + ldb = b_tmp.outerStride(); \ + } \ +\ + BLASPREFIX##hemm_(&side, &uplo, &m, &n, &numext::real_ref(alpha), (const BLASTYPE*)a, &lda, (const BLASTYPE*)b, &ldb, &numext::real_ref(beta), (BLASTYPE*)res, &ldc); \ + } \ +}; + +EIGEN_BLAS_SYMM_R(double, double, d, d) +EIGEN_BLAS_SYMM_R(float, float, f, s) +EIGEN_BLAS_HEMM_R(dcomplex, double, cd, z) +EIGEN_BLAS_HEMM_R(scomplex, float, cf, c) + +} // end namespace internal + +} // end namespace Eigen + +#endif // EIGEN_SELFADJOINT_MATRIX_MATRIX_BLAS_H diff --git a/Eigen/src/Core/products/SelfadjointMatrixMatrix_MKL.h b/Eigen/src/Core/products/SelfadjointMatrixMatrix_MKL.h deleted file mode 100644 index b1176962b..000000000 --- a/Eigen/src/Core/products/SelfadjointMatrixMatrix_MKL.h +++ /dev/null @@ -1,275 +0,0 @@ -/* - 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 Intel(R) MKL - * Self adjoint matrix * matrix product functionality based on ?SYMM/?HEMM. - ******************************************************************************** -*/ - -#ifndef EIGEN_SELFADJOINT_MATRIX_MATRIX_MKL_H -#define EIGEN_SELFADJOINT_MATRIX_MATRIX_MKL_H - -namespace Eigen { - -namespace internal { - - -/* Optimized selfadjoint matrix * matrix (?SYMM/?HEMM) product */ - -#define EIGEN_MKL_SYMM_L(EIGTYPE, BLASTYPE, EIGPREFIX, MKLPREFIX) \ -template \ -struct product_selfadjoint_matrix \ -{\ -\ - static void run( \ - Index rows, Index cols, \ - const EIGTYPE* _lhs, Index lhsStride, \ - const EIGTYPE* _rhs, Index rhsStride, \ - EIGTYPE* res, Index resStride, \ - EIGTYPE alpha, level3_blocking& /*blocking*/) \ - { \ - char side='L', uplo='L'; \ - BlasIndex m, n, lda, ldb, ldc; \ - const EIGTYPE *a, *b; \ - EIGTYPE beta(1); \ - MatrixX##EIGPREFIX b_tmp; \ -\ -/* Set transpose options */ \ -/* Set m, n, k */ \ - m = convert_index(rows); \ - n = convert_index(cols); \ -\ -/* Set lda, ldb, ldc */ \ - lda = convert_index(lhsStride); \ - ldb = convert_index(rhsStride); \ - ldc = convert_index(resStride); \ -\ -/* Set a, b, c */ \ - if (LhsStorageOrder==RowMajor) uplo='U'; \ - a = _lhs; \ -\ - if (RhsStorageOrder==RowMajor) { \ - Map > rhs(_rhs,n,m,OuterStride<>(rhsStride)); \ - b_tmp = rhs.adjoint(); \ - b = b_tmp.data(); \ - ldb = convert_index(b_tmp.outerStride()); \ - } else b = _rhs; \ -\ - MKLPREFIX##symm_(&side, &uplo, &m, &n, &numext::real_ref(alpha), (const BLASTYPE*)a, &lda, (const BLASTYPE*)b, &ldb, &numext::real_ref(beta), (BLASTYPE*)res, &ldc); \ -\ - } \ -}; - - -#define EIGEN_MKL_HEMM_L(EIGTYPE, BLASTYPE, EIGPREFIX, MKLPREFIX) \ -template \ -struct product_selfadjoint_matrix \ -{\ - static void run( \ - Index rows, Index cols, \ - const EIGTYPE* _lhs, Index lhsStride, \ - const EIGTYPE* _rhs, Index rhsStride, \ - EIGTYPE* res, Index resStride, \ - EIGTYPE alpha, level3_blocking& /*blocking*/) \ - { \ - char side='L', uplo='L'; \ - BlasIndex m, n, lda, ldb, ldc; \ - const EIGTYPE *a, *b; \ - EIGTYPE beta(1); \ - MatrixX##EIGPREFIX b_tmp; \ - Matrix a_tmp; \ -\ -/* Set transpose options */ \ -/* Set m, n, k */ \ - m = convert_index(rows); \ - n = convert_index(cols); \ -\ -/* Set lda, ldb, ldc */ \ - lda = convert_index(lhsStride); \ - ldb = convert_index(rhsStride); \ - ldc = convert_index(resStride); \ -\ -/* Set a, b, c */ \ - if (((LhsStorageOrder==ColMajor) && ConjugateLhs) || ((LhsStorageOrder==RowMajor) && (!ConjugateLhs))) { \ - Map, 0, OuterStride<> > lhs(_lhs,m,m,OuterStride<>(lhsStride)); \ - a_tmp = lhs.conjugate(); \ - a = a_tmp.data(); \ - lda = a_tmp.outerStride(); \ - } else a = _lhs; \ - if (LhsStorageOrder==RowMajor) uplo='U'; \ -\ - if (RhsStorageOrder==ColMajor && (!ConjugateRhs)) { \ - b = _rhs; } \ - else { \ - if (RhsStorageOrder==ColMajor && ConjugateRhs) { \ - Map > rhs(_rhs,m,n,OuterStride<>(rhsStride)); \ - b_tmp = rhs.conjugate(); \ - } else \ - if (ConjugateRhs) { \ - Map > rhs(_rhs,n,m,OuterStride<>(rhsStride)); \ - b_tmp = rhs.adjoint(); \ - } else { \ - Map > rhs(_rhs,n,m,OuterStride<>(rhsStride)); \ - b_tmp = rhs.transpose(); \ - } \ - b = b_tmp.data(); \ - ldb = convert_index(b_tmp.outerStride()); \ - } \ -\ - MKLPREFIX##hemm_(&side, &uplo, &m, &n, &numext::real_ref(alpha), (const BLASTYPE*)a, &lda, (const BLASTYPE*)b, &ldb, &numext::real_ref(beta), (BLASTYPE*)res, &ldc); \ -\ - } \ -}; - -EIGEN_MKL_SYMM_L(double, double, d, d) -EIGEN_MKL_SYMM_L(float, float, f, s) -EIGEN_MKL_HEMM_L(dcomplex, double, cd, z) -EIGEN_MKL_HEMM_L(scomplex, float, cf, c) - - -/* Optimized matrix * selfadjoint matrix (?SYMM/?HEMM) product */ - -#define EIGEN_MKL_SYMM_R(EIGTYPE, BLASTYPE, EIGPREFIX, MKLPREFIX) \ -template \ -struct product_selfadjoint_matrix \ -{\ -\ - static void run( \ - Index rows, Index cols, \ - const EIGTYPE* _lhs, Index lhsStride, \ - const EIGTYPE* _rhs, Index rhsStride, \ - EIGTYPE* res, Index resStride, \ - EIGTYPE alpha, level3_blocking& /*blocking*/) \ - { \ - char side='R', uplo='L'; \ - BlasIndex m, n, lda, ldb, ldc; \ - const EIGTYPE *a, *b; \ - EIGTYPE beta(1); \ - MatrixX##EIGPREFIX b_tmp; \ -\ -/* Set m, n, k */ \ - m = convert_index(rows); \ - n = convert_index(cols); \ -\ -/* Set lda, ldb, ldc */ \ - lda = convert_index(rhsStride); \ - ldb = convert_index(lhsStride); \ - ldc = convert_index(resStride); \ -\ -/* Set a, b, c */ \ - if (RhsStorageOrder==RowMajor) uplo='U'; \ - a = _rhs; \ -\ - if (LhsStorageOrder==RowMajor) { \ - Map > lhs(_lhs,n,m,OuterStride<>(rhsStride)); \ - b_tmp = lhs.adjoint(); \ - b = b_tmp.data(); \ - ldb = convert_index(b_tmp.outerStride()); \ - } else b = _lhs; \ -\ - MKLPREFIX##symm_(&side, &uplo, &m, &n, &numext::real_ref(alpha), (const BLASTYPE*)a, &lda, (const BLASTYPE*)b, &ldb, &numext::real_ref(beta), (BLASTYPE*)res, &ldc); \ -\ - } \ -}; - - -#define EIGEN_MKL_HEMM_R(EIGTYPE, BLASTYPE, EIGPREFIX, MKLPREFIX) \ -template \ -struct product_selfadjoint_matrix \ -{\ - static void run( \ - Index rows, Index cols, \ - const EIGTYPE* _lhs, Index lhsStride, \ - const EIGTYPE* _rhs, Index rhsStride, \ - EIGTYPE* res, Index resStride, \ - EIGTYPE alpha, level3_blocking& /*blocking*/) \ - { \ - char side='R', uplo='L'; \ - BlasIndex m, n, lda, ldb, ldc; \ - const EIGTYPE *a, *b; \ - EIGTYPE beta(1); \ - MatrixX##EIGPREFIX b_tmp; \ - Matrix a_tmp; \ -\ -/* Set m, n, k */ \ - m = convert_index(rows); \ - n = convert_index(cols); \ -\ -/* Set lda, ldb, ldc */ \ - lda = convert_index(rhsStride); \ - ldb = convert_index(lhsStride); \ - ldc = convert_index(resStride); \ -\ -/* Set a, b, c */ \ - if (((RhsStorageOrder==ColMajor) && ConjugateRhs) || ((RhsStorageOrder==RowMajor) && (!ConjugateRhs))) { \ - Map, 0, OuterStride<> > rhs(_rhs,n,n,OuterStride<>(rhsStride)); \ - a_tmp = rhs.conjugate(); \ - a = a_tmp.data(); \ - lda = convert_index(a_tmp.outerStride()); \ - } else a = _rhs; \ - if (RhsStorageOrder==RowMajor) uplo='U'; \ -\ - if (LhsStorageOrder==ColMajor && (!ConjugateLhs)) { \ - b = _lhs; } \ - else { \ - if (LhsStorageOrder==ColMajor && ConjugateLhs) { \ - Map > lhs(_lhs,m,n,OuterStride<>(lhsStride)); \ - b_tmp = lhs.conjugate(); \ - } else \ - if (ConjugateLhs) { \ - Map > lhs(_lhs,n,m,OuterStride<>(lhsStride)); \ - b_tmp = lhs.adjoint(); \ - } else { \ - Map > lhs(_lhs,n,m,OuterStride<>(lhsStride)); \ - b_tmp = lhs.transpose(); \ - } \ - b = b_tmp.data(); \ - ldb = b_tmp.outerStride(); \ - } \ -\ - MKLPREFIX##hemm_(&side, &uplo, &m, &n, &numext::real_ref(alpha), (const BLASTYPE*)a, &lda, (const BLASTYPE*)b, &ldb, &numext::real_ref(beta), (BLASTYPE*)res, &ldc); \ - } \ -}; - -EIGEN_MKL_SYMM_R(double, double, d, d) -EIGEN_MKL_SYMM_R(float, float, f, s) -EIGEN_MKL_HEMM_R(dcomplex, double, cd, z) -EIGEN_MKL_HEMM_R(scomplex, float, cf, c) - -} // end namespace internal - -} // end namespace Eigen - -#endif // EIGEN_SELFADJOINT_MATRIX_MATRIX_MKL_H diff --git a/Eigen/src/Core/products/SelfadjointMatrixVector_BLAS.h b/Eigen/src/Core/products/SelfadjointMatrixVector_BLAS.h new file mode 100644 index 000000000..38f23accf --- /dev/null +++ b/Eigen/src/Core/products/SelfadjointMatrixVector_BLAS.h @@ -0,0 +1,111 @@ +/* + 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 + * Selfadjoint matrix-vector product functionality based on ?SYMV/HEMV. + ******************************************************************************** +*/ + +#ifndef EIGEN_SELFADJOINT_MATRIX_VECTOR_BLAS_H +#define EIGEN_SELFADJOINT_MATRIX_VECTOR_BLAS_H + +namespace Eigen { + +namespace internal { + +/********************************************************************** +* This file implements selfadjoint matrix-vector multiplication using BLAS +**********************************************************************/ + +// symv/hemv specialization + +template +struct selfadjoint_matrix_vector_product_symv : + selfadjoint_matrix_vector_product {}; + +#define EIGEN_BLAS_SYMV_SPECIALIZE(Scalar) \ +template \ +struct selfadjoint_matrix_vector_product { \ +static void run( \ + Index size, const Scalar* lhs, Index lhsStride, \ + const Scalar* _rhs, Scalar* res, Scalar alpha) { \ + enum {\ + IsColMajor = StorageOrder==ColMajor \ + }; \ + if (IsColMajor == ConjugateLhs) {\ + selfadjoint_matrix_vector_product::run( \ + size, lhs, lhsStride, _rhs, res, alpha); \ + } else {\ + selfadjoint_matrix_vector_product_symv::run( \ + size, lhs, lhsStride, _rhs, res, alpha); \ + }\ + } \ +}; \ + +EIGEN_BLAS_SYMV_SPECIALIZE(double) +EIGEN_BLAS_SYMV_SPECIALIZE(float) +EIGEN_BLAS_SYMV_SPECIALIZE(dcomplex) +EIGEN_BLAS_SYMV_SPECIALIZE(scomplex) + +#define EIGEN_BLAS_SYMV_SPECIALIZATION(EIGTYPE,BLASTYPE,BLASFUNC) \ +template \ +struct selfadjoint_matrix_vector_product_symv \ +{ \ +typedef Matrix SYMVVector;\ +\ +static void run( \ +Index size, const EIGTYPE* lhs, Index lhsStride, \ +const EIGTYPE* _rhs, EIGTYPE* res, EIGTYPE alpha) \ +{ \ + enum {\ + IsRowMajor = StorageOrder==RowMajor ? 1 : 0, \ + IsLower = UpLo == Lower ? 1 : 0 \ + }; \ + BlasIndex n=convert_index(size), lda=convert_index(lhsStride), incx=1, incy=1; \ + EIGTYPE beta(1); \ + const EIGTYPE *x_ptr; \ + char uplo=(IsRowMajor) ? (IsLower ? 'U' : 'L') : (IsLower ? 'L' : 'U'); \ + SYMVVector x_tmp; \ + if (ConjugateRhs) { \ + Map map_x(_rhs,size,1); \ + x_tmp=map_x.conjugate(); \ + x_ptr=x_tmp.data(); \ + } else x_ptr=_rhs; \ + BLASFUNC(&uplo, &n, &numext::real_ref(alpha), (const BLASTYPE*)lhs, &lda, (const BLASTYPE*)x_ptr, &incx, &numext::real_ref(beta), (BLASTYPE*)res, &incy); \ +}\ +}; + +EIGEN_BLAS_SYMV_SPECIALIZATION(double, double, dsymv_) +EIGEN_BLAS_SYMV_SPECIALIZATION(float, float, ssymv_) +EIGEN_BLAS_SYMV_SPECIALIZATION(dcomplex, double, zhemv_) +EIGEN_BLAS_SYMV_SPECIALIZATION(scomplex, float, chemv_) + +} // end namespace internal + +} // end namespace Eigen + +#endif // EIGEN_SELFADJOINT_MATRIX_VECTOR_BLAS_H diff --git a/Eigen/src/Core/products/SelfadjointMatrixVector_MKL.h b/Eigen/src/Core/products/SelfadjointMatrixVector_MKL.h deleted file mode 100644 index 2a8362202..000000000 --- a/Eigen/src/Core/products/SelfadjointMatrixVector_MKL.h +++ /dev/null @@ -1,111 +0,0 @@ -/* - 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 Intel(R) MKL - * Selfadjoint matrix-vector product functionality based on ?SYMV/HEMV. - ******************************************************************************** -*/ - -#ifndef EIGEN_SELFADJOINT_MATRIX_VECTOR_MKL_H -#define EIGEN_SELFADJOINT_MATRIX_VECTOR_MKL_H - -namespace Eigen { - -namespace internal { - -/********************************************************************** -* This file implements selfadjoint matrix-vector multiplication using BLAS -**********************************************************************/ - -// symv/hemv specialization - -template -struct selfadjoint_matrix_vector_product_symv : - selfadjoint_matrix_vector_product {}; - -#define EIGEN_MKL_SYMV_SPECIALIZE(Scalar) \ -template \ -struct selfadjoint_matrix_vector_product { \ -static void run( \ - Index size, const Scalar* lhs, Index lhsStride, \ - const Scalar* _rhs, Scalar* res, Scalar alpha) { \ - enum {\ - IsColMajor = StorageOrder==ColMajor \ - }; \ - if (IsColMajor == ConjugateLhs) {\ - selfadjoint_matrix_vector_product::run( \ - size, lhs, lhsStride, _rhs, res, alpha); \ - } else {\ - selfadjoint_matrix_vector_product_symv::run( \ - size, lhs, lhsStride, _rhs, res, alpha); \ - }\ - } \ -}; \ - -EIGEN_MKL_SYMV_SPECIALIZE(double) -EIGEN_MKL_SYMV_SPECIALIZE(float) -EIGEN_MKL_SYMV_SPECIALIZE(dcomplex) -EIGEN_MKL_SYMV_SPECIALIZE(scomplex) - -#define EIGEN_MKL_SYMV_SPECIALIZATION(EIGTYPE,BLASTYPE,MKLFUNC) \ -template \ -struct selfadjoint_matrix_vector_product_symv \ -{ \ -typedef Matrix SYMVVector;\ -\ -static void run( \ -Index size, const EIGTYPE* lhs, Index lhsStride, \ -const EIGTYPE* _rhs, EIGTYPE* res, EIGTYPE alpha) \ -{ \ - enum {\ - IsRowMajor = StorageOrder==RowMajor ? 1 : 0, \ - IsLower = UpLo == Lower ? 1 : 0 \ - }; \ - BlasIndex n=convert_index(size), lda=convert_index(lhsStride), incx=1, incy=1; \ - EIGTYPE beta(1); \ - const EIGTYPE *x_ptr; \ - char uplo=(IsRowMajor) ? (IsLower ? 'U' : 'L') : (IsLower ? 'L' : 'U'); \ - SYMVVector x_tmp; \ - if (ConjugateRhs) { \ - Map map_x(_rhs,size,1); \ - x_tmp=map_x.conjugate(); \ - x_ptr=x_tmp.data(); \ - } else x_ptr=_rhs; \ - MKLFUNC(&uplo, &n, &numext::real_ref(alpha), (const BLASTYPE*)lhs, &lda, (const BLASTYPE*)x_ptr, &incx, &numext::real_ref(beta), (BLASTYPE*)res, &incy); \ -}\ -}; - -EIGEN_MKL_SYMV_SPECIALIZATION(double, double, dsymv_) -EIGEN_MKL_SYMV_SPECIALIZATION(float, float, ssymv_) -EIGEN_MKL_SYMV_SPECIALIZATION(dcomplex, double, zhemv_) -EIGEN_MKL_SYMV_SPECIALIZATION(scomplex, float, chemv_) - -} // end namespace internal - -} // end namespace Eigen - -#endif // EIGEN_SELFADJOINT_MATRIX_VECTOR_MKL_H diff --git a/Eigen/src/Core/products/TriangularMatrixMatrix_BLAS.h b/Eigen/src/Core/products/TriangularMatrixMatrix_BLAS.h new file mode 100644 index 000000000..aecded6bb --- /dev/null +++ b/Eigen/src/Core/products/TriangularMatrixMatrix_BLAS.h @@ -0,0 +1,302 @@ +/* + 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_MATRIX_MATRIX_BLAS_H +#define EIGEN_TRIANGULAR_MATRIX_MATRIX_BLAS_H + +namespace Eigen { + +namespace internal { + + +template +struct product_triangular_matrix_matrix_trmm : + product_triangular_matrix_matrix {}; + + +// try to go to BLAS specialization +#define EIGEN_BLAS_TRMM_SPECIALIZE(Scalar, LhsIsTriangular) \ +template \ +struct product_triangular_matrix_matrix { \ + static inline void run(Index _rows, Index _cols, Index _depth, const Scalar* _lhs, Index lhsStride,\ + const Scalar* _rhs, Index rhsStride, Scalar* res, Index resStride, Scalar alpha, level3_blocking& blocking) { \ + product_triangular_matrix_matrix_trmm::run( \ + _rows, _cols, _depth, _lhs, lhsStride, _rhs, rhsStride, res, resStride, alpha, blocking); \ + } \ +}; + +EIGEN_BLAS_TRMM_SPECIALIZE(double, true) +EIGEN_BLAS_TRMM_SPECIALIZE(double, false) +EIGEN_BLAS_TRMM_SPECIALIZE(dcomplex, true) +EIGEN_BLAS_TRMM_SPECIALIZE(dcomplex, false) +EIGEN_BLAS_TRMM_SPECIALIZE(float, true) +EIGEN_BLAS_TRMM_SPECIALIZE(float, false) +EIGEN_BLAS_TRMM_SPECIALIZE(scomplex, true) +EIGEN_BLAS_TRMM_SPECIALIZE(scomplex, false) + +// implements col-major += alpha * op(triangular) * op(general) +#define EIGEN_BLAS_TRMM_L(EIGTYPE, BLASTYPE, EIGPREFIX, BLASPREFIX) \ +template \ +struct product_triangular_matrix_matrix_trmm \ +{ \ + enum { \ + IsLower = (Mode&Lower) == Lower, \ + SetDiag = (Mode&(ZeroDiag|UnitDiag)) ? 0 : 1, \ + IsUnitDiag = (Mode&UnitDiag) ? 1 : 0, \ + IsZeroDiag = (Mode&ZeroDiag) ? 1 : 0, \ + LowUp = IsLower ? Lower : Upper, \ + conjA = ((LhsStorageOrder==ColMajor) && ConjugateLhs) ? 1 : 0 \ + }; \ +\ + static void run( \ + Index _rows, Index _cols, Index _depth, \ + const EIGTYPE* _lhs, Index lhsStride, \ + const EIGTYPE* _rhs, Index rhsStride, \ + EIGTYPE* res, Index resStride, \ + EIGTYPE alpha, level3_blocking& blocking) \ + { \ + Index diagSize = (std::min)(_rows,_depth); \ + Index rows = IsLower ? _rows : diagSize; \ + Index depth = IsLower ? diagSize : _depth; \ + Index cols = _cols; \ +\ + typedef Matrix MatrixLhs; \ + typedef Matrix MatrixRhs; \ +\ +/* Non-square case - doesn't fit to BLAS ?TRMM. Fall to default triangular product or call BLAS ?GEMM*/ \ + if (rows != depth) { \ +\ + /* FIXME handle mkl_domain_get_max_threads */ \ + /*int nthr = mkl_domain_get_max_threads(EIGEN_BLAS_DOMAIN_BLAS);*/ int nthr = 1;\ +\ + if (((nthr==1) && (((std::max)(rows,depth)-diagSize)/(double)diagSize < 0.5))) { \ + /* Most likely no benefit to call TRMM or GEMM from BLAS */ \ + product_triangular_matrix_matrix::run( \ + _rows, _cols, _depth, _lhs, lhsStride, _rhs, rhsStride, res, resStride, alpha, blocking); \ + /*std::cout << "TRMM_L: A is not square! Go to Eigen TRMM implementation!\n";*/ \ + } else { \ + /* Make sense to call GEMM */ \ + Map > lhsMap(_lhs,rows,depth,OuterStride<>(lhsStride)); \ + MatrixLhs aa_tmp=lhsMap.template triangularView(); \ + BlasIndex aStride = convert_index(aa_tmp.outerStride()); \ + gemm_blocking_space gemm_blocking(_rows,_cols,_depth, 1, true); \ + general_matrix_matrix_product::run( \ + rows, cols, depth, aa_tmp.data(), aStride, _rhs, rhsStride, res, resStride, alpha, gemm_blocking, 0); \ +\ + /*std::cout << "TRMM_L: A is not square! Go to BLAS GEMM implementation! " << nthr<<" \n";*/ \ + } \ + return; \ + } \ + char side = 'L', transa, uplo, diag = 'N'; \ + EIGTYPE *b; \ + const EIGTYPE *a; \ + BlasIndex m, n, lda, ldb; \ +\ +/* Set m, n */ \ + m = convert_index(diagSize); \ + n = convert_index(cols); \ +\ +/* Set trans */ \ + transa = (LhsStorageOrder==RowMajor) ? ((ConjugateLhs) ? 'C' : 'T') : 'N'; \ +\ +/* Set b, ldb */ \ + Map > rhs(_rhs,depth,cols,OuterStride<>(rhsStride)); \ + MatrixX##EIGPREFIX b_tmp; \ +\ + if (ConjugateRhs) b_tmp = rhs.conjugate(); else b_tmp = rhs; \ + b = b_tmp.data(); \ + ldb = convert_index(b_tmp.outerStride()); \ +\ +/* Set uplo */ \ + uplo = IsLower ? 'L' : 'U'; \ + if (LhsStorageOrder==RowMajor) uplo = (uplo == 'L') ? 'U' : 'L'; \ +/* Set a, lda */ \ + Map > lhs(_lhs,rows,depth,OuterStride<>(lhsStride)); \ + MatrixLhs a_tmp; \ +\ + if ((conjA!=0) || (SetDiag==0)) { \ + if (conjA) a_tmp = lhs.conjugate(); else a_tmp = lhs; \ + if (IsZeroDiag) \ + a_tmp.diagonal().setZero(); \ + else if (IsUnitDiag) \ + a_tmp.diagonal().setOnes();\ + a = a_tmp.data(); \ + lda = convert_index(a_tmp.outerStride()); \ + } else { \ + a = _lhs; \ + lda = convert_index(lhsStride); \ + } \ + /*std::cout << "TRMM_L: A is square! Go to BLAS TRMM implementation! \n";*/ \ +/* call ?trmm*/ \ + BLASPREFIX##trmm_(&side, &uplo, &transa, &diag, &m, &n, &numext::real_ref(alpha), (const BLASTYPE*)a, &lda, (BLASTYPE*)b, &ldb); \ +\ +/* Add op(a_triangular)*b into res*/ \ + Map > res_tmp(res,rows,cols,OuterStride<>(resStride)); \ + res_tmp=res_tmp+b_tmp; \ + } \ +}; + +EIGEN_BLAS_TRMM_L(double, double, d, d) +EIGEN_BLAS_TRMM_L(dcomplex, double, cd, z) +EIGEN_BLAS_TRMM_L(float, float, f, s) +EIGEN_BLAS_TRMM_L(scomplex, float, cf, c) + +// implements col-major += alpha * op(general) * op(triangular) +#define EIGEN_BLAS_TRMM_R(EIGTYPE, BLASTYPE, EIGPREFIX, BLASPREFIX) \ +template \ +struct product_triangular_matrix_matrix_trmm \ +{ \ + enum { \ + IsLower = (Mode&Lower) == Lower, \ + SetDiag = (Mode&(ZeroDiag|UnitDiag)) ? 0 : 1, \ + IsUnitDiag = (Mode&UnitDiag) ? 1 : 0, \ + IsZeroDiag = (Mode&ZeroDiag) ? 1 : 0, \ + LowUp = IsLower ? Lower : Upper, \ + conjA = ((RhsStorageOrder==ColMajor) && ConjugateRhs) ? 1 : 0 \ + }; \ +\ + static void run( \ + Index _rows, Index _cols, Index _depth, \ + const EIGTYPE* _lhs, Index lhsStride, \ + const EIGTYPE* _rhs, Index rhsStride, \ + EIGTYPE* res, Index resStride, \ + EIGTYPE alpha, level3_blocking& blocking) \ + { \ + Index diagSize = (std::min)(_cols,_depth); \ + Index rows = _rows; \ + Index depth = IsLower ? _depth : diagSize; \ + Index cols = IsLower ? diagSize : _cols; \ +\ + typedef Matrix MatrixLhs; \ + typedef Matrix MatrixRhs; \ +\ +/* Non-square case - doesn't fit to BLAS ?TRMM. Fall to default triangular product or call BLAS ?GEMM*/ \ + if (cols != depth) { \ +\ + int nthr = 1 /*mkl_domain_get_max_threads(EIGEN_BLAS_DOMAIN_BLAS)*/; \ +\ + if ((nthr==1) && (((std::max)(cols,depth)-diagSize)/(double)diagSize < 0.5)) { \ + /* Most likely no benefit to call TRMM or GEMM from BLAS*/ \ + product_triangular_matrix_matrix::run( \ + _rows, _cols, _depth, _lhs, lhsStride, _rhs, rhsStride, res, resStride, alpha, blocking); \ + /*std::cout << "TRMM_R: A is not square! Go to Eigen TRMM implementation!\n";*/ \ + } else { \ + /* Make sense to call GEMM */ \ + Map > rhsMap(_rhs,depth,cols, OuterStride<>(rhsStride)); \ + MatrixRhs aa_tmp=rhsMap.template triangularView(); \ + BlasIndex aStride = convert_index(aa_tmp.outerStride()); \ + gemm_blocking_space gemm_blocking(_rows,_cols,_depth, 1, true); \ + general_matrix_matrix_product::run( \ + rows, cols, depth, _lhs, lhsStride, aa_tmp.data(), aStride, res, resStride, alpha, gemm_blocking, 0); \ +\ + /*std::cout << "TRMM_R: A is not square! Go to BLAS GEMM implementation! " << nthr<<" \n";*/ \ + } \ + return; \ + } \ + char side = 'R', transa, uplo, diag = 'N'; \ + EIGTYPE *b; \ + const EIGTYPE *a; \ + BlasIndex m, n, lda, ldb; \ +\ +/* Set m, n */ \ + m = convert_index(rows); \ + n = convert_index(diagSize); \ +\ +/* Set trans */ \ + transa = (RhsStorageOrder==RowMajor) ? ((ConjugateRhs) ? 'C' : 'T') : 'N'; \ +\ +/* Set b, ldb */ \ + Map > lhs(_lhs,rows,depth,OuterStride<>(lhsStride)); \ + MatrixX##EIGPREFIX b_tmp; \ +\ + if (ConjugateLhs) b_tmp = lhs.conjugate(); else b_tmp = lhs; \ + b = b_tmp.data(); \ + ldb = convert_index(b_tmp.outerStride()); \ +\ +/* Set uplo */ \ + uplo = IsLower ? 'L' : 'U'; \ + if (RhsStorageOrder==RowMajor) uplo = (uplo == 'L') ? 'U' : 'L'; \ +/* Set a, lda */ \ + Map > rhs(_rhs,depth,cols, OuterStride<>(rhsStride)); \ + MatrixRhs a_tmp; \ +\ + if ((conjA!=0) || (SetDiag==0)) { \ + if (conjA) a_tmp = rhs.conjugate(); else a_tmp = rhs; \ + if (IsZeroDiag) \ + a_tmp.diagonal().setZero(); \ + else if (IsUnitDiag) \ + a_tmp.diagonal().setOnes();\ + a = a_tmp.data(); \ + lda = convert_index(a_tmp.outerStride()); \ + } else { \ + a = _rhs; \ + lda = convert_index(rhsStride); \ + } \ + /*std::cout << "TRMM_R: A is square! Go to BLAS TRMM implementation! \n";*/ \ +/* call ?trmm*/ \ + BLASPREFIX##trmm_(&side, &uplo, &transa, &diag, &m, &n, &numext::real_ref(alpha), (const BLASTYPE*)a, &lda, (BLASTYPE*)b, &ldb); \ +\ +/* Add op(a_triangular)*b into res*/ \ + Map > res_tmp(res,rows,cols,OuterStride<>(resStride)); \ + res_tmp=res_tmp+b_tmp; \ + } \ +}; + +EIGEN_BLAS_TRMM_R(double, double, d, d) +EIGEN_BLAS_TRMM_R(dcomplex, double, cd, z) +EIGEN_BLAS_TRMM_R(float, float, f, s) +EIGEN_BLAS_TRMM_R(scomplex, float, cf, c) + +} // end namespace internal + +} // end namespace Eigen + +#endif // EIGEN_TRIANGULAR_MATRIX_MATRIX_BLAS_H diff --git a/Eigen/src/Core/products/TriangularMatrixMatrix_MKL.h b/Eigen/src/Core/products/TriangularMatrixMatrix_MKL.h deleted file mode 100644 index 47a8698a7..000000000 --- a/Eigen/src/Core/products/TriangularMatrixMatrix_MKL.h +++ /dev/null @@ -1,302 +0,0 @@ -/* - 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 Intel(R) MKL - * Triangular matrix * matrix product functionality based on ?TRMM. - ******************************************************************************** -*/ - -#ifndef EIGEN_TRIANGULAR_MATRIX_MATRIX_MKL_H -#define EIGEN_TRIANGULAR_MATRIX_MATRIX_MKL_H - -namespace Eigen { - -namespace internal { - - -template -struct product_triangular_matrix_matrix_trmm : - product_triangular_matrix_matrix {}; - - -// try to go to BLAS specialization -#define EIGEN_MKL_TRMM_SPECIALIZE(Scalar, LhsIsTriangular) \ -template \ -struct product_triangular_matrix_matrix { \ - static inline void run(Index _rows, Index _cols, Index _depth, const Scalar* _lhs, Index lhsStride,\ - const Scalar* _rhs, Index rhsStride, Scalar* res, Index resStride, Scalar alpha, level3_blocking& blocking) { \ - product_triangular_matrix_matrix_trmm::run( \ - _rows, _cols, _depth, _lhs, lhsStride, _rhs, rhsStride, res, resStride, alpha, blocking); \ - } \ -}; - -EIGEN_MKL_TRMM_SPECIALIZE(double, true) -EIGEN_MKL_TRMM_SPECIALIZE(double, false) -EIGEN_MKL_TRMM_SPECIALIZE(dcomplex, true) -EIGEN_MKL_TRMM_SPECIALIZE(dcomplex, false) -EIGEN_MKL_TRMM_SPECIALIZE(float, true) -EIGEN_MKL_TRMM_SPECIALIZE(float, false) -EIGEN_MKL_TRMM_SPECIALIZE(scomplex, true) -EIGEN_MKL_TRMM_SPECIALIZE(scomplex, false) - -// implements col-major += alpha * op(triangular) * op(general) -#define EIGEN_MKL_TRMM_L(EIGTYPE, BLASTYPE, EIGPREFIX, MKLPREFIX) \ -template \ -struct product_triangular_matrix_matrix_trmm \ -{ \ - enum { \ - IsLower = (Mode&Lower) == Lower, \ - SetDiag = (Mode&(ZeroDiag|UnitDiag)) ? 0 : 1, \ - IsUnitDiag = (Mode&UnitDiag) ? 1 : 0, \ - IsZeroDiag = (Mode&ZeroDiag) ? 1 : 0, \ - LowUp = IsLower ? Lower : Upper, \ - conjA = ((LhsStorageOrder==ColMajor) && ConjugateLhs) ? 1 : 0 \ - }; \ -\ - static void run( \ - Index _rows, Index _cols, Index _depth, \ - const EIGTYPE* _lhs, Index lhsStride, \ - const EIGTYPE* _rhs, Index rhsStride, \ - EIGTYPE* res, Index resStride, \ - EIGTYPE alpha, level3_blocking& blocking) \ - { \ - Index diagSize = (std::min)(_rows,_depth); \ - Index rows = IsLower ? _rows : diagSize; \ - Index depth = IsLower ? diagSize : _depth; \ - Index cols = _cols; \ -\ - typedef Matrix MatrixLhs; \ - typedef Matrix MatrixRhs; \ -\ -/* Non-square case - doesn't fit to MKL ?TRMM. Fall to default triangular product or call MKL ?GEMM*/ \ - if (rows != depth) { \ -\ - /* FIXME handle mkl_domain_get_max_threads */ \ - /*int nthr = mkl_domain_get_max_threads(EIGEN_MKL_DOMAIN_BLAS);*/ int nthr = 1;\ -\ - if (((nthr==1) && (((std::max)(rows,depth)-diagSize)/(double)diagSize < 0.5))) { \ - /* Most likely no benefit to call TRMM or GEMM from MKL*/ \ - product_triangular_matrix_matrix::run( \ - _rows, _cols, _depth, _lhs, lhsStride, _rhs, rhsStride, res, resStride, alpha, blocking); \ - /*std::cout << "TRMM_L: A is not square! Go to Eigen TRMM implementation!\n";*/ \ - } else { \ - /* Make sense to call GEMM */ \ - Map > lhsMap(_lhs,rows,depth,OuterStride<>(lhsStride)); \ - MatrixLhs aa_tmp=lhsMap.template triangularView(); \ - BlasIndex aStride = convert_index(aa_tmp.outerStride()); \ - gemm_blocking_space gemm_blocking(_rows,_cols,_depth, 1, true); \ - general_matrix_matrix_product::run( \ - rows, cols, depth, aa_tmp.data(), aStride, _rhs, rhsStride, res, resStride, alpha, gemm_blocking, 0); \ -\ - /*std::cout << "TRMM_L: A is not square! Go to MKL GEMM implementation! " << nthr<<" \n";*/ \ - } \ - return; \ - } \ - char side = 'L', transa, uplo, diag = 'N'; \ - EIGTYPE *b; \ - const EIGTYPE *a; \ - BlasIndex m, n, lda, ldb; \ -\ -/* Set m, n */ \ - m = convert_index(diagSize); \ - n = convert_index(cols); \ -\ -/* Set trans */ \ - transa = (LhsStorageOrder==RowMajor) ? ((ConjugateLhs) ? 'C' : 'T') : 'N'; \ -\ -/* Set b, ldb */ \ - Map > rhs(_rhs,depth,cols,OuterStride<>(rhsStride)); \ - MatrixX##EIGPREFIX b_tmp; \ -\ - if (ConjugateRhs) b_tmp = rhs.conjugate(); else b_tmp = rhs; \ - b = b_tmp.data(); \ - ldb = convert_index(b_tmp.outerStride()); \ -\ -/* Set uplo */ \ - uplo = IsLower ? 'L' : 'U'; \ - if (LhsStorageOrder==RowMajor) uplo = (uplo == 'L') ? 'U' : 'L'; \ -/* Set a, lda */ \ - Map > lhs(_lhs,rows,depth,OuterStride<>(lhsStride)); \ - MatrixLhs a_tmp; \ -\ - if ((conjA!=0) || (SetDiag==0)) { \ - if (conjA) a_tmp = lhs.conjugate(); else a_tmp = lhs; \ - if (IsZeroDiag) \ - a_tmp.diagonal().setZero(); \ - else if (IsUnitDiag) \ - a_tmp.diagonal().setOnes();\ - a = a_tmp.data(); \ - lda = convert_index(a_tmp.outerStride()); \ - } else { \ - a = _lhs; \ - lda = convert_index(lhsStride); \ - } \ - /*std::cout << "TRMM_L: A is square! Go to MKL TRMM implementation! \n";*/ \ -/* call ?trmm*/ \ - MKLPREFIX##trmm_(&side, &uplo, &transa, &diag, &m, &n, &numext::real_ref(alpha), (const BLASTYPE*)a, &lda, (BLASTYPE*)b, &ldb); \ -\ -/* Add op(a_triangular)*b into res*/ \ - Map > res_tmp(res,rows,cols,OuterStride<>(resStride)); \ - res_tmp=res_tmp+b_tmp; \ - } \ -}; - -EIGEN_MKL_TRMM_L(double, double, d, d) -EIGEN_MKL_TRMM_L(dcomplex, double, cd, z) -EIGEN_MKL_TRMM_L(float, float, f, s) -EIGEN_MKL_TRMM_L(scomplex, float, cf, c) - -// implements col-major += alpha * op(general) * op(triangular) -#define EIGEN_MKL_TRMM_R(EIGTYPE, BLASTYPE, EIGPREFIX, MKLPREFIX) \ -template \ -struct product_triangular_matrix_matrix_trmm \ -{ \ - enum { \ - IsLower = (Mode&Lower) == Lower, \ - SetDiag = (Mode&(ZeroDiag|UnitDiag)) ? 0 : 1, \ - IsUnitDiag = (Mode&UnitDiag) ? 1 : 0, \ - IsZeroDiag = (Mode&ZeroDiag) ? 1 : 0, \ - LowUp = IsLower ? Lower : Upper, \ - conjA = ((RhsStorageOrder==ColMajor) && ConjugateRhs) ? 1 : 0 \ - }; \ -\ - static void run( \ - Index _rows, Index _cols, Index _depth, \ - const EIGTYPE* _lhs, Index lhsStride, \ - const EIGTYPE* _rhs, Index rhsStride, \ - EIGTYPE* res, Index resStride, \ - EIGTYPE alpha, level3_blocking& blocking) \ - { \ - Index diagSize = (std::min)(_cols,_depth); \ - Index rows = _rows; \ - Index depth = IsLower ? _depth : diagSize; \ - Index cols = IsLower ? diagSize : _cols; \ -\ - typedef Matrix MatrixLhs; \ - typedef Matrix MatrixRhs; \ -\ -/* Non-square case - doesn't fit to MKL ?TRMM. Fall to default triangular product or call MKL ?GEMM*/ \ - if (cols != depth) { \ -\ - int nthr = 1 /*mkl_domain_get_max_threads(EIGEN_MKL_DOMAIN_BLAS)*/; \ -\ - if ((nthr==1) && (((std::max)(cols,depth)-diagSize)/(double)diagSize < 0.5)) { \ - /* Most likely no benefit to call TRMM or GEMM from MKL*/ \ - product_triangular_matrix_matrix::run( \ - _rows, _cols, _depth, _lhs, lhsStride, _rhs, rhsStride, res, resStride, alpha, blocking); \ - /*std::cout << "TRMM_R: A is not square! Go to Eigen TRMM implementation!\n";*/ \ - } else { \ - /* Make sense to call GEMM */ \ - Map > rhsMap(_rhs,depth,cols, OuterStride<>(rhsStride)); \ - MatrixRhs aa_tmp=rhsMap.template triangularView(); \ - BlasIndex aStride = convert_index(aa_tmp.outerStride()); \ - gemm_blocking_space gemm_blocking(_rows,_cols,_depth, 1, true); \ - general_matrix_matrix_product::run( \ - rows, cols, depth, _lhs, lhsStride, aa_tmp.data(), aStride, res, resStride, alpha, gemm_blocking, 0); \ -\ - /*std::cout << "TRMM_R: A is not square! Go to MKL GEMM implementation! " << nthr<<" \n";*/ \ - } \ - return; \ - } \ - char side = 'R', transa, uplo, diag = 'N'; \ - EIGTYPE *b; \ - const EIGTYPE *a; \ - BlasIndex m, n, lda, ldb; \ -\ -/* Set m, n */ \ - m = convert_index(rows); \ - n = convert_index(diagSize); \ -\ -/* Set trans */ \ - transa = (RhsStorageOrder==RowMajor) ? ((ConjugateRhs) ? 'C' : 'T') : 'N'; \ -\ -/* Set b, ldb */ \ - Map > lhs(_lhs,rows,depth,OuterStride<>(lhsStride)); \ - MatrixX##EIGPREFIX b_tmp; \ -\ - if (ConjugateLhs) b_tmp = lhs.conjugate(); else b_tmp = lhs; \ - b = b_tmp.data(); \ - ldb = convert_index(b_tmp.outerStride()); \ -\ -/* Set uplo */ \ - uplo = IsLower ? 'L' : 'U'; \ - if (RhsStorageOrder==RowMajor) uplo = (uplo == 'L') ? 'U' : 'L'; \ -/* Set a, lda */ \ - Map > rhs(_rhs,depth,cols, OuterStride<>(rhsStride)); \ - MatrixRhs a_tmp; \ -\ - if ((conjA!=0) || (SetDiag==0)) { \ - if (conjA) a_tmp = rhs.conjugate(); else a_tmp = rhs; \ - if (IsZeroDiag) \ - a_tmp.diagonal().setZero(); \ - else if (IsUnitDiag) \ - a_tmp.diagonal().setOnes();\ - a = a_tmp.data(); \ - lda = convert_index(a_tmp.outerStride()); \ - } else { \ - a = _rhs; \ - lda = convert_index(rhsStride); \ - } \ - /*std::cout << "TRMM_R: A is square! Go to MKL TRMM implementation! \n";*/ \ -/* call ?trmm*/ \ - MKLPREFIX##trmm_(&side, &uplo, &transa, &diag, &m, &n, &numext::real_ref(alpha), (const BLASTYPE*)a, &lda, (BLASTYPE*)b, &ldb); \ -\ -/* Add op(a_triangular)*b into res*/ \ - Map > res_tmp(res,rows,cols,OuterStride<>(resStride)); \ - res_tmp=res_tmp+b_tmp; \ - } \ -}; - -EIGEN_MKL_TRMM_R(double, double, d, d) -EIGEN_MKL_TRMM_R(dcomplex, double, cd, z) -EIGEN_MKL_TRMM_R(float, float, f, s) -EIGEN_MKL_TRMM_R(scomplex, float, cf, c) - -} // end namespace internal - -} // end namespace Eigen - -#endif // EIGEN_TRIANGULAR_MATRIX_MATRIX_MKL_H diff --git a/Eigen/src/Core/products/TriangularMatrixVector_BLAS.h b/Eigen/src/Core/products/TriangularMatrixVector_BLAS.h new file mode 100644 index 000000000..07bf26ce5 --- /dev/null +++ b/Eigen/src/Core/products/TriangularMatrixVector_BLAS.h @@ -0,0 +1,241 @@ +/* + 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-vector product functionality based on ?TRMV. + ******************************************************************************** +*/ + +#ifndef EIGEN_TRIANGULAR_MATRIX_VECTOR_BLAS_H +#define EIGEN_TRIANGULAR_MATRIX_VECTOR_BLAS_H + +namespace Eigen { + +namespace internal { + +/********************************************************************** +* This file implements triangular matrix-vector multiplication using BLAS +**********************************************************************/ + +// trmv/hemv specialization + +template +struct triangular_matrix_vector_product_trmv : + triangular_matrix_vector_product {}; + +#define EIGEN_BLAS_TRMV_SPECIALIZE(Scalar) \ +template \ +struct triangular_matrix_vector_product { \ + static void run(Index _rows, Index _cols, const Scalar* _lhs, Index lhsStride, \ + const Scalar* _rhs, Index rhsIncr, Scalar* _res, Index resIncr, Scalar alpha) { \ + triangular_matrix_vector_product_trmv::run( \ + _rows, _cols, _lhs, lhsStride, _rhs, rhsIncr, _res, resIncr, alpha); \ + } \ +}; \ +template \ +struct triangular_matrix_vector_product { \ + static void run(Index _rows, Index _cols, const Scalar* _lhs, Index lhsStride, \ + const Scalar* _rhs, Index rhsIncr, Scalar* _res, Index resIncr, Scalar alpha) { \ + triangular_matrix_vector_product_trmv::run( \ + _rows, _cols, _lhs, lhsStride, _rhs, rhsIncr, _res, resIncr, alpha); \ + } \ +}; + +EIGEN_BLAS_TRMV_SPECIALIZE(double) +EIGEN_BLAS_TRMV_SPECIALIZE(float) +EIGEN_BLAS_TRMV_SPECIALIZE(dcomplex) +EIGEN_BLAS_TRMV_SPECIALIZE(scomplex) + +// implements col-major: res += alpha * op(triangular) * vector +#define EIGEN_BLAS_TRMV_CM(EIGTYPE, BLASTYPE, EIGPREFIX, BLASPREFIX) \ +template \ +struct triangular_matrix_vector_product_trmv { \ + enum { \ + IsLower = (Mode&Lower) == Lower, \ + SetDiag = (Mode&(ZeroDiag|UnitDiag)) ? 0 : 1, \ + IsUnitDiag = (Mode&UnitDiag) ? 1 : 0, \ + IsZeroDiag = (Mode&ZeroDiag) ? 1 : 0, \ + LowUp = IsLower ? Lower : Upper \ + }; \ + static void run(Index _rows, Index _cols, const EIGTYPE* _lhs, Index lhsStride, \ + const EIGTYPE* _rhs, Index rhsIncr, EIGTYPE* _res, Index resIncr, EIGTYPE alpha) \ + { \ + if (ConjLhs || IsZeroDiag) { \ + triangular_matrix_vector_product::run( \ + _rows, _cols, _lhs, lhsStride, _rhs, rhsIncr, _res, resIncr, alpha); \ + return; \ + }\ + Index size = (std::min)(_rows,_cols); \ + Index rows = IsLower ? _rows : size; \ + Index cols = IsLower ? size : _cols; \ +\ + typedef VectorX##EIGPREFIX VectorRhs; \ + EIGTYPE *x, *y;\ +\ +/* Set x*/ \ + Map > rhs(_rhs,cols,InnerStride<>(rhsIncr)); \ + VectorRhs x_tmp; \ + if (ConjRhs) x_tmp = rhs.conjugate(); else x_tmp = rhs; \ + x = x_tmp.data(); \ +\ +/* Square part handling */\ +\ + char trans, uplo, diag; \ + BlasIndex m, n, lda, incx, incy; \ + EIGTYPE const *a; \ + EIGTYPE beta(1); \ +\ +/* Set m, n */ \ + n = convert_index(size); \ + lda = convert_index(lhsStride); \ + incx = 1; \ + incy = convert_index(resIncr); \ +\ +/* Set uplo, trans and diag*/ \ + trans = 'N'; \ + uplo = IsLower ? 'L' : 'U'; \ + diag = IsUnitDiag ? 'U' : 'N'; \ +\ +/* call ?TRMV*/ \ + BLASPREFIX##trmv_(&uplo, &trans, &diag, &n, (const BLASTYPE*)_lhs, &lda, (BLASTYPE*)x, &incx); \ +\ +/* Add op(a_tr)rhs into res*/ \ + BLASPREFIX##axpy_(&n, &numext::real_ref(alpha),(const BLASTYPE*)x, &incx, (BLASTYPE*)_res, &incy); \ +/* Non-square case - doesn't fit to BLAS ?TRMV. Fall to default triangular product*/ \ + if (size<(std::max)(rows,cols)) { \ + if (ConjRhs) x_tmp = rhs.conjugate(); else x_tmp = rhs; \ + x = x_tmp.data(); \ + if (size(rows-size); \ + n = convert_index(size); \ + } \ + else { \ + x += size; \ + y = _res; \ + a = _lhs + size*lda; \ + m = convert_index(size); \ + n = convert_index(cols-size); \ + } \ + BLASPREFIX##gemv_(&trans, &m, &n, &numext::real_ref(alpha), (const BLASTYPE*)a, &lda, (const BLASTYPE*)x, &incx, &numext::real_ref(beta), (BLASTYPE*)y, &incy); \ + } \ + } \ +}; + +EIGEN_BLAS_TRMV_CM(double, double, d, d) +EIGEN_BLAS_TRMV_CM(dcomplex, double, cd, z) +EIGEN_BLAS_TRMV_CM(float, float, f, s) +EIGEN_BLAS_TRMV_CM(scomplex, float, cf, c) + +// implements row-major: res += alpha * op(triangular) * vector +#define EIGEN_BLAS_TRMV_RM(EIGTYPE, BLASTYPE, EIGPREFIX, BLASPREFIX) \ +template \ +struct triangular_matrix_vector_product_trmv { \ + enum { \ + IsLower = (Mode&Lower) == Lower, \ + SetDiag = (Mode&(ZeroDiag|UnitDiag)) ? 0 : 1, \ + IsUnitDiag = (Mode&UnitDiag) ? 1 : 0, \ + IsZeroDiag = (Mode&ZeroDiag) ? 1 : 0, \ + LowUp = IsLower ? Lower : Upper \ + }; \ + static void run(Index _rows, Index _cols, const EIGTYPE* _lhs, Index lhsStride, \ + const EIGTYPE* _rhs, Index rhsIncr, EIGTYPE* _res, Index resIncr, EIGTYPE alpha) \ + { \ + if (IsZeroDiag) { \ + triangular_matrix_vector_product::run( \ + _rows, _cols, _lhs, lhsStride, _rhs, rhsIncr, _res, resIncr, alpha); \ + return; \ + }\ + Index size = (std::min)(_rows,_cols); \ + Index rows = IsLower ? _rows : size; \ + Index cols = IsLower ? size : _cols; \ +\ + typedef VectorX##EIGPREFIX VectorRhs; \ + EIGTYPE *x, *y;\ +\ +/* Set x*/ \ + Map > rhs(_rhs,cols,InnerStride<>(rhsIncr)); \ + VectorRhs x_tmp; \ + if (ConjRhs) x_tmp = rhs.conjugate(); else x_tmp = rhs; \ + x = x_tmp.data(); \ +\ +/* Square part handling */\ +\ + char trans, uplo, diag; \ + BlasIndex m, n, lda, incx, incy; \ + EIGTYPE const *a; \ + EIGTYPE beta(1); \ +\ +/* Set m, n */ \ + n = convert_index(size); \ + lda = convert_index(lhsStride); \ + incx = 1; \ + incy = convert_index(resIncr); \ +\ +/* Set uplo, trans and diag*/ \ + trans = ConjLhs ? 'C' : 'T'; \ + uplo = IsLower ? 'U' : 'L'; \ + diag = IsUnitDiag ? 'U' : 'N'; \ +\ +/* call ?TRMV*/ \ + BLASPREFIX##trmv_(&uplo, &trans, &diag, &n, (const BLASTYPE*)_lhs, &lda, (BLASTYPE*)x, &incx); \ +\ +/* Add op(a_tr)rhs into res*/ \ + BLASPREFIX##axpy_(&n, &numext::real_ref(alpha),(const BLASTYPE*)x, &incx, (BLASTYPE*)_res, &incy); \ +/* Non-square case - doesn't fit to BLAS ?TRMV. Fall to default triangular product*/ \ + if (size<(std::max)(rows,cols)) { \ + if (ConjRhs) x_tmp = rhs.conjugate(); else x_tmp = rhs; \ + x = x_tmp.data(); \ + if (size(rows-size); \ + n = convert_index(size); \ + } \ + else { \ + x += size; \ + y = _res; \ + a = _lhs + size; \ + m = convert_index(size); \ + n = convert_index(cols-size); \ + } \ + BLASPREFIX##gemv_(&trans, &n, &m, &numext::real_ref(alpha), (const BLASTYPE*)a, &lda, (const BLASTYPE*)x, &incx, &numext::real_ref(beta), (BLASTYPE*)y, &incy); \ + } \ + } \ +}; + +EIGEN_BLAS_TRMV_RM(double, double, d, d) +EIGEN_BLAS_TRMV_RM(dcomplex, double, cd, z) +EIGEN_BLAS_TRMV_RM(float, float, f, s) +EIGEN_BLAS_TRMV_RM(scomplex, float, cf, c) + +} // end namespase internal + +} // end namespace Eigen + +#endif // EIGEN_TRIANGULAR_MATRIX_VECTOR_BLAS_H diff --git a/Eigen/src/Core/products/TriangularMatrixVector_MKL.h b/Eigen/src/Core/products/TriangularMatrixVector_MKL.h deleted file mode 100644 index 17c9eeb44..000000000 --- a/Eigen/src/Core/products/TriangularMatrixVector_MKL.h +++ /dev/null @@ -1,241 +0,0 @@ -/* - 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 Intel(R) MKL - * Triangular matrix-vector product functionality based on ?TRMV. - ******************************************************************************** -*/ - -#ifndef EIGEN_TRIANGULAR_MATRIX_VECTOR_MKL_H -#define EIGEN_TRIANGULAR_MATRIX_VECTOR_MKL_H - -namespace Eigen { - -namespace internal { - -/********************************************************************** -* This file implements triangular matrix-vector multiplication using BLAS -**********************************************************************/ - -// trmv/hemv specialization - -template -struct triangular_matrix_vector_product_trmv : - triangular_matrix_vector_product {}; - -#define EIGEN_MKL_TRMV_SPECIALIZE(Scalar) \ -template \ -struct triangular_matrix_vector_product { \ - static void run(Index _rows, Index _cols, const Scalar* _lhs, Index lhsStride, \ - const Scalar* _rhs, Index rhsIncr, Scalar* _res, Index resIncr, Scalar alpha) { \ - triangular_matrix_vector_product_trmv::run( \ - _rows, _cols, _lhs, lhsStride, _rhs, rhsIncr, _res, resIncr, alpha); \ - } \ -}; \ -template \ -struct triangular_matrix_vector_product { \ - static void run(Index _rows, Index _cols, const Scalar* _lhs, Index lhsStride, \ - const Scalar* _rhs, Index rhsIncr, Scalar* _res, Index resIncr, Scalar alpha) { \ - triangular_matrix_vector_product_trmv::run( \ - _rows, _cols, _lhs, lhsStride, _rhs, rhsIncr, _res, resIncr, alpha); \ - } \ -}; - -EIGEN_MKL_TRMV_SPECIALIZE(double) -EIGEN_MKL_TRMV_SPECIALIZE(float) -EIGEN_MKL_TRMV_SPECIALIZE(dcomplex) -EIGEN_MKL_TRMV_SPECIALIZE(scomplex) - -// implements col-major: res += alpha * op(triangular) * vector -#define EIGEN_MKL_TRMV_CM(EIGTYPE, BLASTYPE, EIGPREFIX, MKLPREFIX) \ -template \ -struct triangular_matrix_vector_product_trmv { \ - enum { \ - IsLower = (Mode&Lower) == Lower, \ - SetDiag = (Mode&(ZeroDiag|UnitDiag)) ? 0 : 1, \ - IsUnitDiag = (Mode&UnitDiag) ? 1 : 0, \ - IsZeroDiag = (Mode&ZeroDiag) ? 1 : 0, \ - LowUp = IsLower ? Lower : Upper \ - }; \ - static void run(Index _rows, Index _cols, const EIGTYPE* _lhs, Index lhsStride, \ - const EIGTYPE* _rhs, Index rhsIncr, EIGTYPE* _res, Index resIncr, EIGTYPE alpha) \ - { \ - if (ConjLhs || IsZeroDiag) { \ - triangular_matrix_vector_product::run( \ - _rows, _cols, _lhs, lhsStride, _rhs, rhsIncr, _res, resIncr, alpha); \ - return; \ - }\ - Index size = (std::min)(_rows,_cols); \ - Index rows = IsLower ? _rows : size; \ - Index cols = IsLower ? size : _cols; \ -\ - typedef VectorX##EIGPREFIX VectorRhs; \ - EIGTYPE *x, *y;\ -\ -/* Set x*/ \ - Map > rhs(_rhs,cols,InnerStride<>(rhsIncr)); \ - VectorRhs x_tmp; \ - if (ConjRhs) x_tmp = rhs.conjugate(); else x_tmp = rhs; \ - x = x_tmp.data(); \ -\ -/* Square part handling */\ -\ - char trans, uplo, diag; \ - BlasIndex m, n, lda, incx, incy; \ - EIGTYPE const *a; \ - EIGTYPE beta(1); \ -\ -/* Set m, n */ \ - n = convert_index(size); \ - lda = convert_index(lhsStride); \ - incx = 1; \ - incy = convert_index(resIncr); \ -\ -/* Set uplo, trans and diag*/ \ - trans = 'N'; \ - uplo = IsLower ? 'L' : 'U'; \ - diag = IsUnitDiag ? 'U' : 'N'; \ -\ -/* call ?TRMV*/ \ - MKLPREFIX##trmv_(&uplo, &trans, &diag, &n, (const BLASTYPE*)_lhs, &lda, (BLASTYPE*)x, &incx); \ -\ -/* Add op(a_tr)rhs into res*/ \ - MKLPREFIX##axpy_(&n, &numext::real_ref(alpha),(const BLASTYPE*)x, &incx, (BLASTYPE*)_res, &incy); \ -/* Non-square case - doesn't fit to MKL ?TRMV. Fall to default triangular product*/ \ - if (size<(std::max)(rows,cols)) { \ - if (ConjRhs) x_tmp = rhs.conjugate(); else x_tmp = rhs; \ - x = x_tmp.data(); \ - if (size(rows-size); \ - n = convert_index(size); \ - } \ - else { \ - x += size; \ - y = _res; \ - a = _lhs + size*lda; \ - m = convert_index(size); \ - n = convert_index(cols-size); \ - } \ - MKLPREFIX##gemv_(&trans, &m, &n, &numext::real_ref(alpha), (const BLASTYPE*)a, &lda, (const BLASTYPE*)x, &incx, &numext::real_ref(beta), (BLASTYPE*)y, &incy); \ - } \ - } \ -}; - -EIGEN_MKL_TRMV_CM(double, double, d, d) -EIGEN_MKL_TRMV_CM(dcomplex, double, cd, z) -EIGEN_MKL_TRMV_CM(float, float, f, s) -EIGEN_MKL_TRMV_CM(scomplex, float, cf, c) - -// implements row-major: res += alpha * op(triangular) * vector -#define EIGEN_MKL_TRMV_RM(EIGTYPE, BLASTYPE, EIGPREFIX, MKLPREFIX) \ -template \ -struct triangular_matrix_vector_product_trmv { \ - enum { \ - IsLower = (Mode&Lower) == Lower, \ - SetDiag = (Mode&(ZeroDiag|UnitDiag)) ? 0 : 1, \ - IsUnitDiag = (Mode&UnitDiag) ? 1 : 0, \ - IsZeroDiag = (Mode&ZeroDiag) ? 1 : 0, \ - LowUp = IsLower ? Lower : Upper \ - }; \ - static void run(Index _rows, Index _cols, const EIGTYPE* _lhs, Index lhsStride, \ - const EIGTYPE* _rhs, Index rhsIncr, EIGTYPE* _res, Index resIncr, EIGTYPE alpha) \ - { \ - if (IsZeroDiag) { \ - triangular_matrix_vector_product::run( \ - _rows, _cols, _lhs, lhsStride, _rhs, rhsIncr, _res, resIncr, alpha); \ - return; \ - }\ - Index size = (std::min)(_rows,_cols); \ - Index rows = IsLower ? _rows : size; \ - Index cols = IsLower ? size : _cols; \ -\ - typedef VectorX##EIGPREFIX VectorRhs; \ - EIGTYPE *x, *y;\ -\ -/* Set x*/ \ - Map > rhs(_rhs,cols,InnerStride<>(rhsIncr)); \ - VectorRhs x_tmp; \ - if (ConjRhs) x_tmp = rhs.conjugate(); else x_tmp = rhs; \ - x = x_tmp.data(); \ -\ -/* Square part handling */\ -\ - char trans, uplo, diag; \ - BlasIndex m, n, lda, incx, incy; \ - EIGTYPE const *a; \ - EIGTYPE beta(1); \ -\ -/* Set m, n */ \ - n = convert_index(size); \ - lda = convert_index(lhsStride); \ - incx = 1; \ - incy = convert_index(resIncr); \ -\ -/* Set uplo, trans and diag*/ \ - trans = ConjLhs ? 'C' : 'T'; \ - uplo = IsLower ? 'U' : 'L'; \ - diag = IsUnitDiag ? 'U' : 'N'; \ -\ -/* call ?TRMV*/ \ - MKLPREFIX##trmv_(&uplo, &trans, &diag, &n, (const BLASTYPE*)_lhs, &lda, (BLASTYPE*)x, &incx); \ -\ -/* Add op(a_tr)rhs into res*/ \ - MKLPREFIX##axpy_(&n, &numext::real_ref(alpha),(const BLASTYPE*)x, &incx, (BLASTYPE*)_res, &incy); \ -/* Non-square case - doesn't fit to MKL ?TRMV. Fall to default triangular product*/ \ - if (size<(std::max)(rows,cols)) { \ - if (ConjRhs) x_tmp = rhs.conjugate(); else x_tmp = rhs; \ - x = x_tmp.data(); \ - if (size(rows-size); \ - n = convert_index(size); \ - } \ - else { \ - x += size; \ - y = _res; \ - a = _lhs + size; \ - m = convert_index(size); \ - n = convert_index(cols-size); \ - } \ - MKLPREFIX##gemv_(&trans, &n, &m, &numext::real_ref(alpha), (const BLASTYPE*)a, &lda, (const BLASTYPE*)x, &incx, &numext::real_ref(beta), (BLASTYPE*)y, &incy); \ - } \ - } \ -}; - -EIGEN_MKL_TRMV_RM(double, double, d, d) -EIGEN_MKL_TRMV_RM(dcomplex, double, cd, z) -EIGEN_MKL_TRMV_RM(float, float, f, s) -EIGEN_MKL_TRMV_RM(scomplex, float, cf, c) - -} // end namespase internal - -} // end namespace Eigen - -#endif // EIGEN_TRIANGULAR_MATRIX_VECTOR_MKL_H diff --git a/Eigen/src/Core/products/TriangularSolverMatrix_BLAS.h b/Eigen/src/Core/products/TriangularSolverMatrix_BLAS.h new file mode 100644 index 000000000..88c0fb794 --- /dev/null +++ b/Eigen/src/Core/products/TriangularSolverMatrix_BLAS.h @@ -0,0 +1,151 @@ +/* + 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 + +namespace Eigen { + +namespace internal { + +// implements LeftSide op(triangular)^-1 * general +#define EIGEN_BLAS_TRSM_L(EIGTYPE, BLASTYPE, BLASPREFIX) \ +template \ +struct triangular_solve_matrix \ +{ \ + 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 otherStride, level3_blocking& /*blocking*/) \ + { \ + BlasIndex m = convert_index(size), n = convert_index(otherSize), lda, ldb; \ + char side = 'L', uplo, diag='N', transa; \ + /* Set alpha_ */ \ + EIGTYPE alpha(1); \ + ldb = convert_index(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 MatrixTri; \ + Map > tri(_tri,size,size,OuterStride<>(triStride)); \ + MatrixTri a_tmp; \ +\ + if (conjA) { \ + a_tmp = tri.conjugate(); \ + a = a_tmp.data(); \ + lda = convert_index(a_tmp.outerStride()); \ + } else { \ + a = _tri; \ + lda = convert_index(triStride); \ + } \ + if (IsUnitDiag) diag='U'; \ +/* call ?trsm*/ \ + BLASPREFIX##trsm_(&side, &uplo, &transa, &diag, &m, &n, &numext::real_ref(alpha), (const BLASTYPE*)a, &lda, (BLASTYPE*)_other, &ldb); \ + } \ +}; + +EIGEN_BLAS_TRSM_L(double, double, d) +EIGEN_BLAS_TRSM_L(dcomplex, double, z) +EIGEN_BLAS_TRSM_L(float, float, s) +EIGEN_BLAS_TRSM_L(scomplex, float, c) + + +// implements RightSide general * op(triangular)^-1 +#define EIGEN_BLAS_TRSM_R(EIGTYPE, BLASTYPE, BLASPREFIX) \ +template \ +struct triangular_solve_matrix \ +{ \ + 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 otherStride, level3_blocking& /*blocking*/) \ + { \ + BlasIndex m = convert_index(otherSize), n = convert_index(size), lda, ldb; \ + char side = 'R', uplo, diag='N', transa; \ + /* Set alpha_ */ \ + EIGTYPE alpha(1); \ + ldb = convert_index(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 MatrixTri; \ + Map > tri(_tri,size,size,OuterStride<>(triStride)); \ + MatrixTri a_tmp; \ +\ + if (conjA) { \ + a_tmp = tri.conjugate(); \ + a = a_tmp.data(); \ + lda = convert_index(a_tmp.outerStride()); \ + } else { \ + a = _tri; \ + lda = convert_index(triStride); \ + } \ + if (IsUnitDiag) diag='U'; \ +/* call ?trsm*/ \ + BLASPREFIX##trsm_(&side, &uplo, &transa, &diag, &m, &n, &numext::real_ref(alpha), (const BLASTYPE*)a, &lda, (BLASTYPE*)_other, &ldb); \ + /*std::cout << "TRMS_L specialization!\n";*/ \ + } \ +}; + +EIGEN_BLAS_TRSM_R(double, double, d) +EIGEN_BLAS_TRSM_R(dcomplex, double, z) +EIGEN_BLAS_TRSM_R(float, float, s) +EIGEN_BLAS_TRSM_R(scomplex, float, c) + + +} // end namespace internal + +} // end namespace Eigen + +#endif // EIGEN_TRIANGULAR_SOLVER_MATRIX_BLAS_H diff --git a/Eigen/src/Core/products/TriangularSolverMatrix_MKL.h b/Eigen/src/Core/products/TriangularSolverMatrix_MKL.h deleted file mode 100644 index 1f68a1cec..000000000 --- a/Eigen/src/Core/products/TriangularSolverMatrix_MKL.h +++ /dev/null @@ -1,151 +0,0 @@ -/* - 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 Intel(R) MKL - * Triangular matrix * matrix product functionality based on ?TRMM. - ******************************************************************************** -*/ - -#ifndef EIGEN_TRIANGULAR_SOLVER_MATRIX_MKL_H -#define EIGEN_TRIANGULAR_SOLVER_MATRIX_MKL_H - -namespace Eigen { - -namespace internal { - -// implements LeftSide op(triangular)^-1 * general -#define EIGEN_MKL_TRSM_L(EIGTYPE, BLASTYPE, MKLPREFIX) \ -template \ -struct triangular_solve_matrix \ -{ \ - 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 otherStride, level3_blocking& /*blocking*/) \ - { \ - BlasIndex m = convert_index(size), n = convert_index(otherSize), lda, ldb; \ - char side = 'L', uplo, diag='N', transa; \ - /* Set alpha_ */ \ - EIGTYPE alpha(1); \ - ldb = convert_index(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 MatrixTri; \ - Map > tri(_tri,size,size,OuterStride<>(triStride)); \ - MatrixTri a_tmp; \ -\ - if (conjA) { \ - a_tmp = tri.conjugate(); \ - a = a_tmp.data(); \ - lda = convert_index(a_tmp.outerStride()); \ - } else { \ - a = _tri; \ - lda = convert_index(triStride); \ - } \ - if (IsUnitDiag) diag='U'; \ -/* call ?trsm*/ \ - MKLPREFIX##trsm_(&side, &uplo, &transa, &diag, &m, &n, &numext::real_ref(alpha), (const BLASTYPE*)a, &lda, (BLASTYPE*)_other, &ldb); \ - } \ -}; - -EIGEN_MKL_TRSM_L(double, double, d) -EIGEN_MKL_TRSM_L(dcomplex, double, z) -EIGEN_MKL_TRSM_L(float, float, s) -EIGEN_MKL_TRSM_L(scomplex, float, c) - - -// implements RightSide general * op(triangular)^-1 -#define EIGEN_MKL_TRSM_R(EIGTYPE, BLASTYPE, MKLPREFIX) \ -template \ -struct triangular_solve_matrix \ -{ \ - 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 otherStride, level3_blocking& /*blocking*/) \ - { \ - BlasIndex m = convert_index(otherSize), n = convert_index(size), lda, ldb; \ - char side = 'R', uplo, diag='N', transa; \ - /* Set alpha_ */ \ - EIGTYPE alpha(1); \ - ldb = convert_index(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 MatrixTri; \ - Map > tri(_tri,size,size,OuterStride<>(triStride)); \ - MatrixTri a_tmp; \ -\ - if (conjA) { \ - a_tmp = tri.conjugate(); \ - a = a_tmp.data(); \ - lda = convert_index(a_tmp.outerStride()); \ - } else { \ - a = _tri; \ - lda = convert_index(triStride); \ - } \ - if (IsUnitDiag) diag='U'; \ -/* call ?trsm*/ \ - MKLPREFIX##trsm_(&side, &uplo, &transa, &diag, &m, &n, &numext::real_ref(alpha), (const BLASTYPE*)a, &lda, (BLASTYPE*)_other, &ldb); \ - /*std::cout << "TRMS_L specialization!\n";*/ \ - } \ -}; - -EIGEN_MKL_TRSM_R(double, double, d) -EIGEN_MKL_TRSM_R(dcomplex, double, z) -EIGEN_MKL_TRSM_R(float, float, s) -EIGEN_MKL_TRSM_R(scomplex, float, c) - - -} // end namespace internal - -} // end namespace Eigen - -#endif // EIGEN_TRIANGULAR_SOLVER_MATRIX_MKL_H -- cgit v1.2.3