acado: SelfadjointMatrixMatrix

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  *   Self adjoint matrix * matrix product functionality based on ?SYMM/?HEMM.
00030  ********************************************************************************
00031 */
00032 
00033 #ifndef EIGEN_SELFADJOINT_MATRIX_MATRIX_MKL_H
00034 #define EIGEN_SELFADJOINT_MATRIX_MATRIX_MKL_H
00035 
00036 namespace Eigen { 
00037 
00038 namespace internal {
00039 
00040 
00041 /* Optimized selfadjoint matrix * matrix (?SYMM/?HEMM) product */
00042 
00043 #define EIGEN_MKL_SYMM_L(EIGTYPE, MKLTYPE, EIGPREFIX, MKLPREFIX) \
00044 template <typename Index, \
00045           int LhsStorageOrder, bool ConjugateLhs, \
00046           int RhsStorageOrder, bool ConjugateRhs> \
00047 struct product_selfadjoint_matrix<EIGTYPE,Index,LhsStorageOrder,true,ConjugateLhs,RhsStorageOrder,false,ConjugateRhs,ColMajor> \
00048 {\
00049 \
00050   static void run( \
00051     Index rows, Index cols, \
00052     const EIGTYPE* _lhs, Index lhsStride, \
00053     const EIGTYPE* _rhs, Index rhsStride, \
00054     EIGTYPE* res,        Index resStride, \
00055     EIGTYPE alpha) \
00056   { \
00057     char side='L', uplo='L'; \
00058     MKL_INT m, n, lda, ldb, ldc; \
00059     const EIGTYPE *a, *b; \
00060     MKLTYPE alpha_, beta_; \
00061     MatrixX##EIGPREFIX b_tmp; \
00062     EIGTYPE myone(1);\
00063 \
00064 /* Set transpose options */ \
00065 /* Set m, n, k */ \
00066     m = (MKL_INT)rows;  \
00067     n = (MKL_INT)cols;  \
00068 \
00069 /* Set alpha_ & beta_ */ \
00070     assign_scalar_eig2mkl(alpha_, alpha); \
00071     assign_scalar_eig2mkl(beta_, myone); \
00072 \
00073 /* Set lda, ldb, ldc */ \
00074     lda = (MKL_INT)lhsStride; \
00075     ldb = (MKL_INT)rhsStride; \
00076     ldc = (MKL_INT)resStride; \
00077 \
00078 /* Set a, b, c */ \
00079     if (LhsStorageOrder==RowMajor) uplo='U'; \
00080     a = _lhs; \
00081 \
00082     if (RhsStorageOrder==RowMajor) { \
00083       Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > rhs(_rhs,n,m,OuterStride<>(rhsStride)); \
00084       b_tmp = rhs.adjoint(); \
00085       b = b_tmp.data(); \
00086       ldb = b_tmp.outerStride(); \
00087     } else b = _rhs; \
00088 \
00089     MKLPREFIX##symm(&side, &uplo, &m, &n, &alpha_, (const MKLTYPE*)a, &lda, (const MKLTYPE*)b, &ldb, &beta_, (MKLTYPE*)res, &ldc); \
00090 \
00091   } \
00092 };
00093 
00094 
00095 #define EIGEN_MKL_HEMM_L(EIGTYPE, MKLTYPE, EIGPREFIX, MKLPREFIX) \
00096 template <typename Index, \
00097           int LhsStorageOrder, bool ConjugateLhs, \
00098           int RhsStorageOrder, bool ConjugateRhs> \
00099 struct product_selfadjoint_matrix<EIGTYPE,Index,LhsStorageOrder,true,ConjugateLhs,RhsStorageOrder,false,ConjugateRhs,ColMajor> \
00100 {\
00101   static void run( \
00102     Index rows, Index cols, \
00103     const EIGTYPE* _lhs, Index lhsStride, \
00104     const EIGTYPE* _rhs, Index rhsStride, \
00105     EIGTYPE* res,        Index resStride, \
00106     EIGTYPE alpha) \
00107   { \
00108     char side='L', uplo='L'; \
00109     MKL_INT m, n, lda, ldb, ldc; \
00110     const EIGTYPE *a, *b; \
00111     MKLTYPE alpha_, beta_; \
00112     MatrixX##EIGPREFIX b_tmp; \
00113     Matrix<EIGTYPE, Dynamic, Dynamic, LhsStorageOrder> a_tmp; \
00114     EIGTYPE myone(1); \
00115 \
00116 /* Set transpose options */ \
00117 /* Set m, n, k */ \
00118     m = (MKL_INT)rows; \
00119     n = (MKL_INT)cols; \
00120 \
00121 /* Set alpha_ & beta_ */ \
00122     assign_scalar_eig2mkl(alpha_, alpha); \
00123     assign_scalar_eig2mkl(beta_, myone); \
00124 \
00125 /* Set lda, ldb, ldc */ \
00126     lda = (MKL_INT)lhsStride; \
00127     ldb = (MKL_INT)rhsStride; \
00128     ldc = (MKL_INT)resStride; \
00129 \
00130 /* Set a, b, c */ \
00131     if (((LhsStorageOrder==ColMajor) && ConjugateLhs) || ((LhsStorageOrder==RowMajor) && (!ConjugateLhs))) { \
00132       Map<const Matrix<EIGTYPE, Dynamic, Dynamic, LhsStorageOrder>, 0, OuterStride<> > lhs(_lhs,m,m,OuterStride<>(lhsStride)); \
00133       a_tmp = lhs.conjugate(); \
00134       a = a_tmp.data(); \
00135       lda = a_tmp.outerStride(); \
00136     } else a = _lhs; \
00137     if (LhsStorageOrder==RowMajor) uplo='U'; \
00138 \
00139     if (RhsStorageOrder==ColMajor && (!ConjugateRhs)) { \
00140        b = _rhs; } \
00141     else { \
00142       if (RhsStorageOrder==ColMajor && ConjugateRhs) { \
00143         Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > rhs(_rhs,m,n,OuterStride<>(rhsStride)); \
00144         b_tmp = rhs.conjugate(); \
00145       } else \
00146       if (ConjugateRhs) { \
00147         Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > rhs(_rhs,n,m,OuterStride<>(rhsStride)); \
00148         b_tmp = rhs.adjoint(); \
00149       } else { \
00150         Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > rhs(_rhs,n,m,OuterStride<>(rhsStride)); \
00151         b_tmp = rhs.transpose(); \
00152       } \
00153       b = b_tmp.data(); \
00154       ldb = b_tmp.outerStride(); \
00155     } \
00156 \
00157     MKLPREFIX##hemm(&side, &uplo, &m, &n, &alpha_, (const MKLTYPE*)a, &lda, (const MKLTYPE*)b, &ldb, &beta_, (MKLTYPE*)res, &ldc); \
00158 \
00159   } \
00160 };
00161 
00162 EIGEN_MKL_SYMM_L(double, double, d, d)
00163 EIGEN_MKL_SYMM_L(float, float, f, s)
00164 EIGEN_MKL_HEMM_L(dcomplex, MKL_Complex16, cd, z)
00165 EIGEN_MKL_HEMM_L(scomplex, MKL_Complex8, cf, c)
00166 
00167 
00168 /* Optimized matrix * selfadjoint matrix (?SYMM/?HEMM) product */
00169 
00170 #define EIGEN_MKL_SYMM_R(EIGTYPE, MKLTYPE, EIGPREFIX, MKLPREFIX) \
00171 template <typename Index, \
00172           int LhsStorageOrder, bool ConjugateLhs, \
00173           int RhsStorageOrder, bool ConjugateRhs> \
00174 struct product_selfadjoint_matrix<EIGTYPE,Index,LhsStorageOrder,false,ConjugateLhs,RhsStorageOrder,true,ConjugateRhs,ColMajor> \
00175 {\
00176 \
00177   static void run( \
00178     Index rows, Index cols, \
00179     const EIGTYPE* _lhs, Index lhsStride, \
00180     const EIGTYPE* _rhs, Index rhsStride, \
00181     EIGTYPE* res,        Index resStride, \
00182     EIGTYPE alpha) \
00183   { \
00184     char side='R', uplo='L'; \
00185     MKL_INT m, n, lda, ldb, ldc; \
00186     const EIGTYPE *a, *b; \
00187     MKLTYPE alpha_, beta_; \
00188     MatrixX##EIGPREFIX b_tmp; \
00189     EIGTYPE myone(1);\
00190 \
00191 /* Set m, n, k */ \
00192     m = (MKL_INT)rows;  \
00193     n = (MKL_INT)cols;  \
00194 \
00195 /* Set alpha_ & beta_ */ \
00196     assign_scalar_eig2mkl(alpha_, alpha); \
00197     assign_scalar_eig2mkl(beta_, myone); \
00198 \
00199 /* Set lda, ldb, ldc */ \
00200     lda = (MKL_INT)rhsStride; \
00201     ldb = (MKL_INT)lhsStride; \
00202     ldc = (MKL_INT)resStride; \
00203 \
00204 /* Set a, b, c */ \
00205     if (RhsStorageOrder==RowMajor) uplo='U'; \
00206     a = _rhs; \
00207 \
00208     if (LhsStorageOrder==RowMajor) { \
00209       Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > lhs(_lhs,n,m,OuterStride<>(rhsStride)); \
00210       b_tmp = lhs.adjoint(); \
00211       b = b_tmp.data(); \
00212       ldb = b_tmp.outerStride(); \
00213     } else b = _lhs; \
00214 \
00215     MKLPREFIX##symm(&side, &uplo, &m, &n, &alpha_, (const MKLTYPE*)a, &lda, (const MKLTYPE*)b, &ldb, &beta_, (MKLTYPE*)res, &ldc); \
00216 \
00217   } \
00218 };
00219 
00220 
00221 #define EIGEN_MKL_HEMM_R(EIGTYPE, MKLTYPE, EIGPREFIX, MKLPREFIX) \
00222 template <typename Index, \
00223           int LhsStorageOrder, bool ConjugateLhs, \
00224           int RhsStorageOrder, bool ConjugateRhs> \
00225 struct product_selfadjoint_matrix<EIGTYPE,Index,LhsStorageOrder,false,ConjugateLhs,RhsStorageOrder,true,ConjugateRhs,ColMajor> \
00226 {\
00227   static void run( \
00228     Index rows, Index cols, \
00229     const EIGTYPE* _lhs, Index lhsStride, \
00230     const EIGTYPE* _rhs, Index rhsStride, \
00231     EIGTYPE* res,        Index resStride, \
00232     EIGTYPE alpha) \
00233   { \
00234     char side='R', uplo='L'; \
00235     MKL_INT m, n, lda, ldb, ldc; \
00236     const EIGTYPE *a, *b; \
00237     MKLTYPE alpha_, beta_; \
00238     MatrixX##EIGPREFIX b_tmp; \
00239     Matrix<EIGTYPE, Dynamic, Dynamic, RhsStorageOrder> a_tmp; \
00240     EIGTYPE myone(1); \
00241 \
00242 /* Set m, n, k */ \
00243     m = (MKL_INT)rows; \
00244     n = (MKL_INT)cols; \
00245 \
00246 /* Set alpha_ & beta_ */ \
00247     assign_scalar_eig2mkl(alpha_, alpha); \
00248     assign_scalar_eig2mkl(beta_, myone); \
00249 \
00250 /* Set lda, ldb, ldc */ \
00251     lda = (MKL_INT)rhsStride; \
00252     ldb = (MKL_INT)lhsStride; \
00253     ldc = (MKL_INT)resStride; \
00254 \
00255 /* Set a, b, c */ \
00256     if (((RhsStorageOrder==ColMajor) && ConjugateRhs) || ((RhsStorageOrder==RowMajor) && (!ConjugateRhs))) { \
00257       Map<const Matrix<EIGTYPE, Dynamic, Dynamic, RhsStorageOrder>, 0, OuterStride<> > rhs(_rhs,n,n,OuterStride<>(rhsStride)); \
00258       a_tmp = rhs.conjugate(); \
00259       a = a_tmp.data(); \
00260       lda = a_tmp.outerStride(); \
00261     } else a = _rhs; \
00262     if (RhsStorageOrder==RowMajor) uplo='U'; \
00263 \
00264     if (LhsStorageOrder==ColMajor && (!ConjugateLhs)) { \
00265        b = _lhs; } \
00266     else { \
00267       if (LhsStorageOrder==ColMajor && ConjugateLhs) { \
00268         Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > lhs(_lhs,m,n,OuterStride<>(lhsStride)); \
00269         b_tmp = lhs.conjugate(); \
00270       } else \
00271       if (ConjugateLhs) { \
00272         Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > lhs(_lhs,n,m,OuterStride<>(lhsStride)); \
00273         b_tmp = lhs.adjoint(); \
00274       } else { \
00275         Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > lhs(_lhs,n,m,OuterStride<>(lhsStride)); \
00276         b_tmp = lhs.transpose(); \
00277       } \
00278       b = b_tmp.data(); \
00279       ldb = b_tmp.outerStride(); \
00280     } \
00281 \
00282     MKLPREFIX##hemm(&side, &uplo, &m, &n, &alpha_, (const MKLTYPE*)a, &lda, (const MKLTYPE*)b, &ldb, &beta_, (MKLTYPE*)res, &ldc); \
00283   } \
00284 };
00285 
00286 EIGEN_MKL_SYMM_R(double, double, d, d)
00287 EIGEN_MKL_SYMM_R(float, float, f, s)
00288 EIGEN_MKL_HEMM_R(dcomplex, MKL_Complex16, cd, z)
00289 EIGEN_MKL_HEMM_R(scomplex, MKL_Complex8, cf, c)
00290 
00291 } // end namespace internal
00292 
00293 } // end namespace Eigen
00294 
00295 #endif // EIGEN_SELFADJOINT_MATRIX_MATRIX_MKL_H