Eigne Matrix Class Library
Dependents: Eigen_test Odometry_test AttitudeEstimation_usingTicker MPU9250_Quaternion_Binary_Serial ... more
SelfadjointMatrixMatrix_MKL.h
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
Generated on Tue Jul 12 2022 17:47:00 by 1.7.2