TriangularSolverMatrix_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  *   Triangular matrix * matrix product functionality based on ?TRMM.
00030  ********************************************************************************
00031 */
00032 
00033 #ifndef EIGEN_TRIANGULAR_SOLVER_MATRIX_MKL_H
00034 #define EIGEN_TRIANGULAR_SOLVER_MATRIX_MKL_H
00035 
00036 namespace Eigen {
00037 
00038 namespace internal {
00039 
00040 // implements LeftSide op(triangular)^-1 * general
00041 #define EIGEN_MKL_TRSM_L(EIGTYPE, MKLTYPE, MKLPREFIX) \
00042 template <typename Index, int Mode, bool Conjugate, int TriStorageOrder> \
00043 struct triangular_solve_matrix<EIGTYPE,Index,OnTheLeft,Mode,Conjugate,TriStorageOrder,ColMajor> \
00044 { \
00045   enum { \
00046     IsLower = (Mode&Lower) == Lower, \
00047     IsUnitDiag  = (Mode&UnitDiag) ? 1 : 0, \
00048     IsZeroDiag  = (Mode&ZeroDiag) ? 1 : 0, \
00049     conjA = ((TriStorageOrder==ColMajor) && Conjugate) ? 1 : 0 \
00050   }; \
00051   static void run( \
00052       Index size, Index otherSize, \
00053       const EIGTYPE* _tri, Index triStride, \
00054       EIGTYPE* _other, Index otherStride, level3_blocking<EIGTYPE,EIGTYPE>& /*blocking*/) \
00055   { \
00056    MKL_INT m = size, n = otherSize, lda, ldb; \
00057    char side = 'L', uplo, diag='N', transa; \
00058    /* Set alpha_ */ \
00059    MKLTYPE alpha; \
00060    EIGTYPE myone(1); \
00061    assign_scalar_eig2mkl(alpha, myone); \
00062    ldb = otherStride;\
00063 \
00064    const EIGTYPE *a; \
00065 /* Set trans */ \
00066    transa = (TriStorageOrder==RowMajor) ? ((Conjugate) ? 'C' : 'T') : 'N'; \
00067 /* Set uplo */ \
00068    uplo = IsLower ? 'L' : 'U'; \
00069    if (TriStorageOrder==RowMajor) uplo = (uplo == 'L') ? 'U' : 'L'; \
00070 /* Set a, lda */ \
00071    typedef Matrix<EIGTYPE, Dynamic, Dynamic, TriStorageOrder> MatrixTri; \
00072    Map<const MatrixTri, 0, OuterStride<> > tri(_tri,size,size,OuterStride<>(triStride)); \
00073    MatrixTri a_tmp; \
00074 \
00075    if (conjA) { \
00076      a_tmp = tri.conjugate(); \
00077      a = a_tmp.data(); \
00078      lda = a_tmp.outerStride(); \
00079    } else { \
00080      a = _tri; \
00081      lda = triStride; \
00082    } \
00083    if (IsUnitDiag) diag='U'; \
00084 /* call ?trsm*/ \
00085    MKLPREFIX##trsm(&side, &uplo, &transa, &diag, &m, &n, &alpha, (const MKLTYPE*)a, &lda, (MKLTYPE*)_other, &ldb); \
00086  } \
00087 };
00088 
00089 EIGEN_MKL_TRSM_L(double, double, d)
00090 EIGEN_MKL_TRSM_L(dcomplex, MKL_Complex16, z)
00091 EIGEN_MKL_TRSM_L(float, float, s)
00092 EIGEN_MKL_TRSM_L(scomplex, MKL_Complex8, c)
00093 
00094 
00095 // implements RightSide general * op(triangular)^-1
00096 #define EIGEN_MKL_TRSM_R(EIGTYPE, MKLTYPE, MKLPREFIX) \
00097 template <typename Index, int Mode, bool Conjugate, int TriStorageOrder> \
00098 struct triangular_solve_matrix<EIGTYPE,Index,OnTheRight,Mode,Conjugate,TriStorageOrder,ColMajor> \
00099 { \
00100   enum { \
00101     IsLower = (Mode&Lower) == Lower, \
00102     IsUnitDiag  = (Mode&UnitDiag) ? 1 : 0, \
00103     IsZeroDiag  = (Mode&ZeroDiag) ? 1 : 0, \
00104     conjA = ((TriStorageOrder==ColMajor) && Conjugate) ? 1 : 0 \
00105   }; \
00106   static void run( \
00107       Index size, Index otherSize, \
00108       const EIGTYPE* _tri, Index triStride, \
00109       EIGTYPE* _other, Index otherStride, level3_blocking<EIGTYPE,EIGTYPE>& /*blocking*/) \
00110   { \
00111    MKL_INT m = otherSize, n = size, lda, ldb; \
00112    char side = 'R', uplo, diag='N', transa; \
00113    /* Set alpha_ */ \
00114    MKLTYPE alpha; \
00115    EIGTYPE myone(1); \
00116    assign_scalar_eig2mkl(alpha, myone); \
00117    ldb = otherStride;\
00118 \
00119    const EIGTYPE *a; \
00120 /* Set trans */ \
00121    transa = (TriStorageOrder==RowMajor) ? ((Conjugate) ? 'C' : 'T') : 'N'; \
00122 /* Set uplo */ \
00123    uplo = IsLower ? 'L' : 'U'; \
00124    if (TriStorageOrder==RowMajor) uplo = (uplo == 'L') ? 'U' : 'L'; \
00125 /* Set a, lda */ \
00126    typedef Matrix<EIGTYPE, Dynamic, Dynamic, TriStorageOrder> MatrixTri; \
00127    Map<const MatrixTri, 0, OuterStride<> > tri(_tri,size,size,OuterStride<>(triStride)); \
00128    MatrixTri a_tmp; \
00129 \
00130    if (conjA) { \
00131      a_tmp = tri.conjugate(); \
00132      a = a_tmp.data(); \
00133      lda = a_tmp.outerStride(); \
00134    } else { \
00135      a = _tri; \
00136      lda = triStride; \
00137    } \
00138    if (IsUnitDiag) diag='U'; \
00139 /* call ?trsm*/ \
00140    MKLPREFIX##trsm(&side, &uplo, &transa, &diag, &m, &n, &alpha, (const MKLTYPE*)a, &lda, (MKLTYPE*)_other, &ldb); \
00141    /*std::cout << "TRMS_L specialization!\n";*/ \
00142  } \
00143 };
00144 
00145 EIGEN_MKL_TRSM_R(double, double, d)
00146 EIGEN_MKL_TRSM_R(dcomplex, MKL_Complex16, z)
00147 EIGEN_MKL_TRSM_R(float, float, s)
00148 EIGEN_MKL_TRSM_R(scomplex, MKL_Complex8, c)
00149 
00150 
00151 } // end namespace internal
00152 
00153 } // end namespace Eigen
00154 
00155 #endif // EIGEN_TRIANGULAR_SOLVER_MATRIX_MKL_H


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Thu Aug 27 2015 12:01:17