1 /* 2 Copyright (c) 2011, Intel Corporation. All rights reserved. 3 4 Redistribution and use in source and binary forms, with or without modification, 5 are permitted provided that the following conditions are met: 6 7 * Redistributions of source code must retain the above copyright notice, this 8 list of conditions and the following disclaimer. 9 * Redistributions in binary form must reproduce the above copyright notice, 10 this list of conditions and the following disclaimer in the documentation 11 and/or other materials provided with the distribution. 12 * Neither the name of Intel Corporation nor the names of its contributors may 13 be used to endorse or promote products derived from this software without 14 specific prior written permission. 15 16 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND 17 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 18 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 19 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR 20 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES 21 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 22 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON 23 ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 24 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 25 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 26 27 ******************************************************************************** 28 * Content : Eigen bindings to Intel(R) MKL 29 * Self adjoint matrix * matrix product functionality based on ?SYMM/?HEMM. 30 ******************************************************************************** 31 */ 32 33 #ifndef EIGEN_SELFADJOINT_MATRIX_MATRIX_MKL_H 34 #define EIGEN_SELFADJOINT_MATRIX_MATRIX_MKL_H 35 36 namespace Eigen { 37 38 namespace internal { 39 40 41 /* Optimized selfadjoint matrix * matrix (?SYMM/?HEMM) product */ 42 43 #define EIGEN_MKL_SYMM_L(EIGTYPE, MKLTYPE, EIGPREFIX, MKLPREFIX) \ 44 template <typename Index, \ 45 int LhsStorageOrder, bool ConjugateLhs, \ 46 int RhsStorageOrder, bool ConjugateRhs> \ 47 struct product_selfadjoint_matrix<EIGTYPE,Index,LhsStorageOrder,true,ConjugateLhs,RhsStorageOrder,false,ConjugateRhs,ColMajor> \ 48 {\ 49 \ 50 static EIGEN_DONT_INLINE void run( \ 51 Index rows, Index cols, \ 52 const EIGTYPE* _lhs, Index lhsStride, \ 53 const EIGTYPE* _rhs, Index rhsStride, \ 54 EIGTYPE* res, Index resStride, \ 55 EIGTYPE alpha) \ 56 { \ 57 char side='L', uplo='L'; \ 58 MKL_INT m, n, lda, ldb, ldc; \ 59 const EIGTYPE *a, *b; \ 60 MKLTYPE alpha_, beta_; \ 61 MatrixX##EIGPREFIX b_tmp; \ 62 EIGTYPE myone(1);\ 63 \ 64 /* Set transpose options */ \ 65 /* Set m, n, k */ \ 66 m = (MKL_INT)rows; \ 67 n = (MKL_INT)cols; \ 68 \ 69 /* Set alpha_ & beta_ */ \ 70 assign_scalar_eig2mkl(alpha_, alpha); \ 71 assign_scalar_eig2mkl(beta_, myone); \ 72 \ 73 /* Set lda, ldb, ldc */ \ 74 lda = (MKL_INT)lhsStride; \ 75 ldb = (MKL_INT)rhsStride; \ 76 ldc = (MKL_INT)resStride; \ 77 \ 78 /* Set a, b, c */ \ 79 if (LhsStorageOrder==RowMajor) uplo='U'; \ 80 a = _lhs; \ 81 \ 82 if (RhsStorageOrder==RowMajor) { \ 83 Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > rhs(_rhs,n,m,OuterStride<>(rhsStride)); \ 84 b_tmp = rhs.adjoint(); \ 85 b = b_tmp.data(); \ 86 ldb = b_tmp.outerStride(); \ 87 } else b = _rhs; \ 88 \ 89 MKLPREFIX##symm(&side, &uplo, &m, &n, &alpha_, (const MKLTYPE*)a, &lda, (const MKLTYPE*)b, &ldb, &beta_, (MKLTYPE*)res, &ldc); \ 90 \ 91 } \ 92 }; 93 94 95 #define EIGEN_MKL_HEMM_L(EIGTYPE, MKLTYPE, EIGPREFIX, MKLPREFIX) \ 96 template <typename Index, \ 97 int LhsStorageOrder, bool ConjugateLhs, \ 98 int RhsStorageOrder, bool ConjugateRhs> \ 99 struct product_selfadjoint_matrix<EIGTYPE,Index,LhsStorageOrder,true,ConjugateLhs,RhsStorageOrder,false,ConjugateRhs,ColMajor> \ 100 {\ 101 static EIGEN_DONT_INLINE void run( \ 102 Index rows, Index cols, \ 103 const EIGTYPE* _lhs, Index lhsStride, \ 104 const EIGTYPE* _rhs, Index rhsStride, \ 105 EIGTYPE* res, Index resStride, \ 106 EIGTYPE alpha) \ 107 { \ 108 char side='L', uplo='L'; \ 109 MKL_INT m, n, lda, ldb, ldc; \ 110 const EIGTYPE *a, *b; \ 111 MKLTYPE alpha_, beta_; \ 112 MatrixX##EIGPREFIX b_tmp; \ 113 Matrix<EIGTYPE, Dynamic, Dynamic, LhsStorageOrder> a_tmp; \ 114 EIGTYPE myone(1); \ 115 \ 116 /* Set transpose options */ \ 117 /* Set m, n, k */ \ 118 m = (MKL_INT)rows; \ 119 n = (MKL_INT)cols; \ 120 \ 121 /* Set alpha_ & beta_ */ \ 122 assign_scalar_eig2mkl(alpha_, alpha); \ 123 assign_scalar_eig2mkl(beta_, myone); \ 124 \ 125 /* Set lda, ldb, ldc */ \ 126 lda = (MKL_INT)lhsStride; \ 127 ldb = (MKL_INT)rhsStride; \ 128 ldc = (MKL_INT)resStride; \ 129 \ 130 /* Set a, b, c */ \ 131 if (((LhsStorageOrder==ColMajor) && ConjugateLhs) || ((LhsStorageOrder==RowMajor) && (!ConjugateLhs))) { \ 132 Map<const Matrix<EIGTYPE, Dynamic, Dynamic, LhsStorageOrder>, 0, OuterStride<> > lhs(_lhs,m,m,OuterStride<>(lhsStride)); \ 133 a_tmp = lhs.conjugate(); \ 134 a = a_tmp.data(); \ 135 lda = a_tmp.outerStride(); \ 136 } else a = _lhs; \ 137 if (LhsStorageOrder==RowMajor) uplo='U'; \ 138 \ 139 if (RhsStorageOrder==ColMajor && (!ConjugateRhs)) { \ 140 b = _rhs; } \ 141 else { \ 142 if (RhsStorageOrder==ColMajor && ConjugateRhs) { \ 143 Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > rhs(_rhs,m,n,OuterStride<>(rhsStride)); \ 144 b_tmp = rhs.conjugate(); \ 145 } else \ 146 if (ConjugateRhs) { \ 147 Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > rhs(_rhs,n,m,OuterStride<>(rhsStride)); \ 148 b_tmp = rhs.adjoint(); \ 149 } else { \ 150 Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > rhs(_rhs,n,m,OuterStride<>(rhsStride)); \ 151 b_tmp = rhs.transpose(); \ 152 } \ 153 b = b_tmp.data(); \ 154 ldb = b_tmp.outerStride(); \ 155 } \ 156 \ 157 MKLPREFIX##hemm(&side, &uplo, &m, &n, &alpha_, (const MKLTYPE*)a, &lda, (const MKLTYPE*)b, &ldb, &beta_, (MKLTYPE*)res, &ldc); \ 158 \ 159 } \ 160 }; 161 162 EIGEN_MKL_SYMM_L(double, double, d, d) 163 EIGEN_MKL_SYMM_L(float, float, f, s) 164 EIGEN_MKL_HEMM_L(dcomplex, MKL_Complex16, cd, z) 165 EIGEN_MKL_HEMM_L(scomplex, MKL_Complex8, cf, c) 166 167 168 /* Optimized matrix * selfadjoint matrix (?SYMM/?HEMM) product */ 169 170 #define EIGEN_MKL_SYMM_R(EIGTYPE, MKLTYPE, EIGPREFIX, MKLPREFIX) \ 171 template <typename Index, \ 172 int LhsStorageOrder, bool ConjugateLhs, \ 173 int RhsStorageOrder, bool ConjugateRhs> \ 174 struct product_selfadjoint_matrix<EIGTYPE,Index,LhsStorageOrder,false,ConjugateLhs,RhsStorageOrder,true,ConjugateRhs,ColMajor> \ 175 {\ 176 \ 177 static EIGEN_DONT_INLINE void run( \ 178 Index rows, Index cols, \ 179 const EIGTYPE* _lhs, Index lhsStride, \ 180 const EIGTYPE* _rhs, Index rhsStride, \ 181 EIGTYPE* res, Index resStride, \ 182 EIGTYPE alpha) \ 183 { \ 184 char side='R', uplo='L'; \ 185 MKL_INT m, n, lda, ldb, ldc; \ 186 const EIGTYPE *a, *b; \ 187 MKLTYPE alpha_, beta_; \ 188 MatrixX##EIGPREFIX b_tmp; \ 189 EIGTYPE myone(1);\ 190 \ 191 /* Set m, n, k */ \ 192 m = (MKL_INT)rows; \ 193 n = (MKL_INT)cols; \ 194 \ 195 /* Set alpha_ & beta_ */ \ 196 assign_scalar_eig2mkl(alpha_, alpha); \ 197 assign_scalar_eig2mkl(beta_, myone); \ 198 \ 199 /* Set lda, ldb, ldc */ \ 200 lda = (MKL_INT)rhsStride; \ 201 ldb = (MKL_INT)lhsStride; \ 202 ldc = (MKL_INT)resStride; \ 203 \ 204 /* Set a, b, c */ \ 205 if (RhsStorageOrder==RowMajor) uplo='U'; \ 206 a = _rhs; \ 207 \ 208 if (LhsStorageOrder==RowMajor) { \ 209 Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > lhs(_lhs,n,m,OuterStride<>(rhsStride)); \ 210 b_tmp = lhs.adjoint(); \ 211 b = b_tmp.data(); \ 212 ldb = b_tmp.outerStride(); \ 213 } else b = _lhs; \ 214 \ 215 MKLPREFIX##symm(&side, &uplo, &m, &n, &alpha_, (const MKLTYPE*)a, &lda, (const MKLTYPE*)b, &ldb, &beta_, (MKLTYPE*)res, &ldc); \ 216 \ 217 } \ 218 }; 219 220 221 #define EIGEN_MKL_HEMM_R(EIGTYPE, MKLTYPE, EIGPREFIX, MKLPREFIX) \ 222 template <typename Index, \ 223 int LhsStorageOrder, bool ConjugateLhs, \ 224 int RhsStorageOrder, bool ConjugateRhs> \ 225 struct product_selfadjoint_matrix<EIGTYPE,Index,LhsStorageOrder,false,ConjugateLhs,RhsStorageOrder,true,ConjugateRhs,ColMajor> \ 226 {\ 227 static EIGEN_DONT_INLINE void run( \ 228 Index rows, Index cols, \ 229 const EIGTYPE* _lhs, Index lhsStride, \ 230 const EIGTYPE* _rhs, Index rhsStride, \ 231 EIGTYPE* res, Index resStride, \ 232 EIGTYPE alpha) \ 233 { \ 234 char side='R', uplo='L'; \ 235 MKL_INT m, n, lda, ldb, ldc; \ 236 const EIGTYPE *a, *b; \ 237 MKLTYPE alpha_, beta_; \ 238 MatrixX##EIGPREFIX b_tmp; \ 239 Matrix<EIGTYPE, Dynamic, Dynamic, RhsStorageOrder> a_tmp; \ 240 EIGTYPE myone(1); \ 241 \ 242 /* Set m, n, k */ \ 243 m = (MKL_INT)rows; \ 244 n = (MKL_INT)cols; \ 245 \ 246 /* Set alpha_ & beta_ */ \ 247 assign_scalar_eig2mkl(alpha_, alpha); \ 248 assign_scalar_eig2mkl(beta_, myone); \ 249 \ 250 /* Set lda, ldb, ldc */ \ 251 lda = (MKL_INT)rhsStride; \ 252 ldb = (MKL_INT)lhsStride; \ 253 ldc = (MKL_INT)resStride; \ 254 \ 255 /* Set a, b, c */ \ 256 if (((RhsStorageOrder==ColMajor) && ConjugateRhs) || ((RhsStorageOrder==RowMajor) && (!ConjugateRhs))) { \ 257 Map<const Matrix<EIGTYPE, Dynamic, Dynamic, RhsStorageOrder>, 0, OuterStride<> > rhs(_rhs,n,n,OuterStride<>(rhsStride)); \ 258 a_tmp = rhs.conjugate(); \ 259 a = a_tmp.data(); \ 260 lda = a_tmp.outerStride(); \ 261 } else a = _rhs; \ 262 if (RhsStorageOrder==RowMajor) uplo='U'; \ 263 \ 264 if (LhsStorageOrder==ColMajor && (!ConjugateLhs)) { \ 265 b = _lhs; } \ 266 else { \ 267 if (LhsStorageOrder==ColMajor && ConjugateLhs) { \ 268 Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > lhs(_lhs,m,n,OuterStride<>(lhsStride)); \ 269 b_tmp = lhs.conjugate(); \ 270 } else \ 271 if (ConjugateLhs) { \ 272 Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > lhs(_lhs,n,m,OuterStride<>(lhsStride)); \ 273 b_tmp = lhs.adjoint(); \ 274 } else { \ 275 Map<const MatrixX##EIGPREFIX, 0, OuterStride<> > lhs(_lhs,n,m,OuterStride<>(lhsStride)); \ 276 b_tmp = lhs.transpose(); \ 277 } \ 278 b = b_tmp.data(); \ 279 ldb = b_tmp.outerStride(); \ 280 } \ 281 \ 282 MKLPREFIX##hemm(&side, &uplo, &m, &n, &alpha_, (const MKLTYPE*)a, &lda, (const MKLTYPE*)b, &ldb, &beta_, (MKLTYPE*)res, &ldc); \ 283 } \ 284 }; 285 286 EIGEN_MKL_SYMM_R(double, double, d, d) 287 EIGEN_MKL_SYMM_R(float, float, f, s) 288 EIGEN_MKL_HEMM_R(dcomplex, MKL_Complex16, cd, z) 289 EIGEN_MKL_HEMM_R(scomplex, MKL_Complex8, cf, c) 290 291 } // end namespace internal 292 293 } // end namespace Eigen 294 295 #endif // EIGEN_SELFADJOINT_MATRIX_MATRIX_MKL_H 296
