GeneralMatrixMatrixTriangular_MKL.h
Go to the documentation of this file.
00001 /*
00002  Copyright (c) 2011, Intel Corporation. All rights reserved.
00003 
00004  Redistribution and use in source and binary forms, with or without modification,
00005  are permitted provided that the following conditions are met:
00006 
00007  * Redistributions of source code must retain the above copyright notice, this
00008    list of conditions and the following disclaimer.
00009  * Redistributions in binary form must reproduce the above copyright notice,
00010    this list of conditions and the following disclaimer in the documentation
00011    and/or other materials provided with the distribution.
00012  * Neither the name of Intel Corporation nor the names of its contributors may
00013    be used to endorse or promote products derived from this software without
00014    specific prior written permission.
00015 
00016  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
00017  ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
00018  WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
00019  DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
00020  ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
00021  (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
00022  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
00023  ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
00024  (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00025  SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00026 
00027  ********************************************************************************
00028  *   Content : Eigen bindings to Intel(R) MKL
00029  *   Level 3 BLAS SYRK/HERK implementation.
00030  ********************************************************************************
00031 */
00032 
00033 #ifndef EIGEN_GENERAL_MATRIX_MATRIX_TRIANGULAR_MKL_H
00034 #define EIGEN_GENERAL_MATRIX_MATRIX_TRIANGULAR_MKL_H
00035 
00036 namespace Eigen { 
00037 
00038 namespace internal {
00039 
00040 template <typename Index, typename Scalar, int AStorageOrder, bool ConjugateA, int ResStorageOrder, int  UpLo>
00041 struct general_matrix_matrix_rankupdate :
00042        general_matrix_matrix_triangular_product<
00043          Index,Scalar,AStorageOrder,ConjugateA,Scalar,AStorageOrder,ConjugateA,ResStorageOrder,UpLo,BuiltIn> {};
00044 
00045 
00046 // try to go to BLAS specialization
00047 #define EIGEN_MKL_RANKUPDATE_SPECIALIZE(Scalar) \
00048 template <typename Index, int LhsStorageOrder, bool ConjugateLhs, \
00049                           int RhsStorageOrder, bool ConjugateRhs, int  UpLo> \
00050 struct general_matrix_matrix_triangular_product<Index,Scalar,LhsStorageOrder,ConjugateLhs, \
00051                Scalar,RhsStorageOrder,ConjugateRhs,ColMajor,UpLo,Specialized> { \
00052   static EIGEN_STRONG_INLINE void run(Index size, Index depth,const Scalar* lhs, Index lhsStride, \
00053                           const Scalar* rhs, Index rhsStride, Scalar* res, Index resStride, Scalar alpha) \
00054   { \
00055     if (lhs==rhs) { \
00056       general_matrix_matrix_rankupdate<Index,Scalar,LhsStorageOrder,ConjugateLhs,ColMajor,UpLo> \
00057       ::run(size,depth,lhs,lhsStride,rhs,rhsStride,res,resStride,alpha); \
00058     } else { \
00059       general_matrix_matrix_triangular_product<Index, \
00060         Scalar, LhsStorageOrder, ConjugateLhs, \
00061         Scalar, RhsStorageOrder, ConjugateRhs, \
00062         ColMajor, UpLo, BuiltIn> \
00063       ::run(size,depth,lhs,lhsStride,rhs,rhsStride,res,resStride,alpha); \
00064     } \
00065   } \
00066 };
00067 
00068 EIGEN_MKL_RANKUPDATE_SPECIALIZE(double)
00069 //EIGEN_MKL_RANKUPDATE_SPECIALIZE(dcomplex)
00070 EIGEN_MKL_RANKUPDATE_SPECIALIZE(float)
00071 //EIGEN_MKL_RANKUPDATE_SPECIALIZE(scomplex)
00072 
00073 // SYRK for float/double
00074 #define EIGEN_MKL_RANKUPDATE_R(EIGTYPE, MKLTYPE, MKLFUNC) \
00075 template <typename Index, int AStorageOrder, bool ConjugateA, int  UpLo> \
00076 struct general_matrix_matrix_rankupdate<Index,EIGTYPE,AStorageOrder,ConjugateA,ColMajor,UpLo> { \
00077   enum { \
00078     IsLower = (UpLo&Lower) == Lower, \
00079     LowUp = IsLower ? Lower : Upper, \
00080     conjA = ((AStorageOrder==ColMajor) && ConjugateA) ? 1 : 0 \
00081   }; \
00082   static EIGEN_STRONG_INLINE void run(Index size, Index depth,const EIGTYPE* lhs, Index lhsStride, \
00083                           const EIGTYPE* rhs, Index rhsStride, EIGTYPE* res, Index resStride, EIGTYPE alpha) \
00084   { \
00085   /* typedef Matrix<EIGTYPE, Dynamic, Dynamic, RhsStorageOrder> MatrixRhs;*/ \
00086 \
00087    MKL_INT lda=lhsStride, ldc=resStride, n=size, k=depth; \
00088    char uplo=(IsLower) ? 'L' : 'U', trans=(AStorageOrder==RowMajor) ? 'T':'N'; \
00089    MKLTYPE alpha_, beta_; \
00090 \
00091 /* Set alpha_ & beta_ */ \
00092    assign_scalar_eig2mkl<MKLTYPE, EIGTYPE>(alpha_, alpha); \
00093    assign_scalar_eig2mkl<MKLTYPE, EIGTYPE>(beta_, EIGTYPE(1)); \
00094    MKLFUNC(&uplo, &trans, &n, &k, &alpha_, lhs, &lda, &beta_, res, &ldc); \
00095   } \
00096 };
00097 
00098 // HERK for complex data
00099 #define EIGEN_MKL_RANKUPDATE_C(EIGTYPE, MKLTYPE, RTYPE, MKLFUNC) \
00100 template <typename Index, int AStorageOrder, bool ConjugateA, int  UpLo> \
00101 struct general_matrix_matrix_rankupdate<Index,EIGTYPE,AStorageOrder,ConjugateA,ColMajor,UpLo> { \
00102   enum { \
00103     IsLower = (UpLo&Lower) == Lower, \
00104     LowUp = IsLower ? Lower : Upper, \
00105     conjA = (((AStorageOrder==ColMajor) && ConjugateA) || ((AStorageOrder==RowMajor) && !ConjugateA)) ? 1 : 0 \
00106   }; \
00107   static EIGEN_STRONG_INLINE void run(Index size, Index depth,const EIGTYPE* lhs, Index lhsStride, \
00108                           const EIGTYPE* rhs, Index rhsStride, EIGTYPE* res, Index resStride, EIGTYPE alpha) \
00109   { \
00110    typedef Matrix<EIGTYPE, Dynamic, Dynamic, AStorageOrder> MatrixType; \
00111 \
00112    MKL_INT lda=lhsStride, ldc=resStride, n=size, k=depth; \
00113    char uplo=(IsLower) ? 'L' : 'U', trans=(AStorageOrder==RowMajor) ? 'C':'N'; \
00114    RTYPE alpha_, beta_; \
00115    const EIGTYPE* a_ptr; \
00116 \
00117 /* Set alpha_ & beta_ */ \
00118 /*   assign_scalar_eig2mkl<MKLTYPE, EIGTYPE>(alpha_, alpha); */\
00119 /*   assign_scalar_eig2mkl<MKLTYPE, EIGTYPE>(beta_, EIGTYPE(1));*/ \
00120    alpha_ = alpha.real(); \
00121    beta_ = 1.0; \
00122 /* Copy with conjugation in some cases*/ \
00123    MatrixType a; \
00124    if (conjA) { \
00125      Map<const MatrixType, 0, OuterStride<> > mapA(lhs,n,k,OuterStride<>(lhsStride)); \
00126      a = mapA.conjugate(); \
00127      lda = a.outerStride(); \
00128      a_ptr = a.data(); \
00129    } else a_ptr=lhs; \
00130    MKLFUNC(&uplo, &trans, &n, &k, &alpha_, (MKLTYPE*)a_ptr, &lda, &beta_, (MKLTYPE*)res, &ldc); \
00131   } \
00132 };
00133 
00134 
00135 EIGEN_MKL_RANKUPDATE_R(double, double, dsyrk)
00136 EIGEN_MKL_RANKUPDATE_R(float,  float,  ssyrk)
00137 
00138 //EIGEN_MKL_RANKUPDATE_C(dcomplex, MKL_Complex16, double, zherk)
00139 //EIGEN_MKL_RANKUPDATE_C(scomplex, MKL_Complex8,  double, cherk)
00140 
00141 
00142 } // end namespace internal
00143 
00144 } // end namespace Eigen
00145 
00146 #endif // EIGEN_GENERAL_MATRIX_MATRIX_TRIANGULAR_MKL_H


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:37:14