level3_impl.h revision 7faaa9f3f0df9d23790277834d426c3d992ac3ba
1// This file is part of Eigen, a lightweight C++ template library 2// for linear algebra. 3// 4// Copyright (C) 2009-2010 Gael Guennebaud <gael.guennebaud@inria.fr> 5// 6// This Source Code Form is subject to the terms of the Mozilla 7// Public License v. 2.0. If a copy of the MPL was not distributed 8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10#include "common.h" 11 12int EIGEN_BLAS_FUNC(gemm)(char *opa, char *opb, int *m, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc) 13{ 14// std::cerr << "in gemm " << *opa << " " << *opb << " " << *m << " " << *n << " " << *k << " " << *lda << " " << *ldb << " " << *ldc << " " << *palpha << " " << *pbeta << "\n"; 15 typedef void (*functype)(DenseIndex, DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, Scalar, internal::level3_blocking<Scalar,Scalar>&, Eigen::internal::GemmParallelInfo<DenseIndex>*); 16 static functype func[12]; 17 18 static bool init = false; 19 if(!init) 20 { 21 for(int k=0; k<12; ++k) 22 func[k] = 0; 23 func[NOTR | (NOTR << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,ColMajor,false,Scalar,ColMajor,false,ColMajor>::run); 24 func[TR | (NOTR << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,false,Scalar,ColMajor,false,ColMajor>::run); 25 func[ADJ | (NOTR << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,false,ColMajor>::run); 26 func[NOTR | (TR << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,false,ColMajor>::run); 27 func[TR | (TR << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,false,Scalar,RowMajor,false,ColMajor>::run); 28 func[ADJ | (TR << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,RowMajor,false,ColMajor>::run); 29 func[NOTR | (ADJ << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,Conj, ColMajor>::run); 30 func[TR | (ADJ << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,false,Scalar,RowMajor,Conj, ColMajor>::run); 31 func[ADJ | (ADJ << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,RowMajor,Conj, ColMajor>::run); 32 init = true; 33 } 34 35 Scalar* a = reinterpret_cast<Scalar*>(pa); 36 Scalar* b = reinterpret_cast<Scalar*>(pb); 37 Scalar* c = reinterpret_cast<Scalar*>(pc); 38 Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 39 Scalar beta = *reinterpret_cast<Scalar*>(pbeta); 40 41 int info = 0; 42 if(OP(*opa)==INVALID) info = 1; 43 else if(OP(*opb)==INVALID) info = 2; 44 else if(*m<0) info = 3; 45 else if(*n<0) info = 4; 46 else if(*k<0) info = 5; 47 else if(*lda<std::max(1,(OP(*opa)==NOTR)?*m:*k)) info = 8; 48 else if(*ldb<std::max(1,(OP(*opb)==NOTR)?*k:*n)) info = 10; 49 else if(*ldc<std::max(1,*m)) info = 13; 50 if(info) 51 return xerbla_(SCALAR_SUFFIX_UP"GEMM ",&info,6); 52 53 if(beta!=Scalar(1)) 54 { 55 if(beta==Scalar(0)) matrix(c, *m, *n, *ldc).setZero(); 56 else matrix(c, *m, *n, *ldc) *= beta; 57 } 58 59 internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic> blocking(*m,*n,*k); 60 61 int code = OP(*opa) | (OP(*opb) << 2); 62 func[code](*m, *n, *k, a, *lda, b, *ldb, c, *ldc, alpha, blocking, 0); 63 return 0; 64} 65 66int EIGEN_BLAS_FUNC(trsm)(char *side, char *uplo, char *opa, char *diag, int *m, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb) 67{ 68// std::cerr << "in trsm " << *side << " " << *uplo << " " << *opa << " " << *diag << " " << *m << "," << *n << " " << *palpha << " " << *lda << " " << *ldb<< "\n"; 69 typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, internal::level3_blocking<Scalar,Scalar>&); 70 static functype func[32]; 71 72 static bool init = false; 73 if(!init) 74 { 75 for(int k=0; k<32; ++k) 76 func[k] = 0; 77 78 func[NOTR | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|0, false,ColMajor,ColMajor>::run); 79 func[TR | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|0, false,RowMajor,ColMajor>::run); 80 func[ADJ | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|0, Conj, RowMajor,ColMajor>::run); 81 82 func[NOTR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|0, false,ColMajor,ColMajor>::run); 83 func[TR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|0, false,RowMajor,ColMajor>::run); 84 func[ADJ | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|0, Conj, RowMajor,ColMajor>::run); 85 86 func[NOTR | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|0, false,ColMajor,ColMajor>::run); 87 func[TR | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|0, false,RowMajor,ColMajor>::run); 88 func[ADJ | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|0, Conj, RowMajor,ColMajor>::run); 89 90 func[NOTR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|0, false,ColMajor,ColMajor>::run); 91 func[TR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|0, false,RowMajor,ColMajor>::run); 92 func[ADJ | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|0, Conj, RowMajor,ColMajor>::run); 93 94 95 func[NOTR | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|UnitDiag,false,ColMajor,ColMajor>::run); 96 func[TR | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|UnitDiag,false,RowMajor,ColMajor>::run); 97 func[ADJ | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|UnitDiag,Conj, RowMajor,ColMajor>::run); 98 99 func[NOTR | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|UnitDiag,false,ColMajor,ColMajor>::run); 100 func[TR | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|UnitDiag,false,RowMajor,ColMajor>::run); 101 func[ADJ | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|UnitDiag,Conj, RowMajor,ColMajor>::run); 102 103 func[NOTR | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|UnitDiag,false,ColMajor,ColMajor>::run); 104 func[TR | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|UnitDiag,false,RowMajor,ColMajor>::run); 105 func[ADJ | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|UnitDiag,Conj, RowMajor,ColMajor>::run); 106 107 func[NOTR | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|UnitDiag,false,ColMajor,ColMajor>::run); 108 func[TR | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|UnitDiag,false,RowMajor,ColMajor>::run); 109 func[ADJ | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|UnitDiag,Conj, RowMajor,ColMajor>::run); 110 111 init = true; 112 } 113 114 Scalar* a = reinterpret_cast<Scalar*>(pa); 115 Scalar* b = reinterpret_cast<Scalar*>(pb); 116 Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 117 118 int info = 0; 119 if(SIDE(*side)==INVALID) info = 1; 120 else if(UPLO(*uplo)==INVALID) info = 2; 121 else if(OP(*opa)==INVALID) info = 3; 122 else if(DIAG(*diag)==INVALID) info = 4; 123 else if(*m<0) info = 5; 124 else if(*n<0) info = 6; 125 else if(*lda<std::max(1,(SIDE(*side)==LEFT)?*m:*n)) info = 9; 126 else if(*ldb<std::max(1,*m)) info = 11; 127 if(info) 128 return xerbla_(SCALAR_SUFFIX_UP"TRSM ",&info,6); 129 130 int code = OP(*opa) | (SIDE(*side) << 2) | (UPLO(*uplo) << 3) | (DIAG(*diag) << 4); 131 132 if(SIDE(*side)==LEFT) 133 { 134 internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic,4> blocking(*m,*n,*m); 135 func[code](*m, *n, a, *lda, b, *ldb, blocking); 136 } 137 else 138 { 139 internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic,4> blocking(*m,*n,*n); 140 func[code](*n, *m, a, *lda, b, *ldb, blocking); 141 } 142 143 if(alpha!=Scalar(1)) 144 matrix(b,*m,*n,*ldb) *= alpha; 145 146 return 0; 147} 148 149 150// b = alpha*op(a)*b for side = 'L'or'l' 151// b = alpha*b*op(a) for side = 'R'or'r' 152int EIGEN_BLAS_FUNC(trmm)(char *side, char *uplo, char *opa, char *diag, int *m, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb) 153{ 154// std::cerr << "in trmm " << *side << " " << *uplo << " " << *opa << " " << *diag << " " << *m << " " << *n << " " << *lda << " " << *ldb << " " << *palpha << "\n"; 155 typedef void (*functype)(DenseIndex, DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, const Scalar&, internal::level3_blocking<Scalar,Scalar>&); 156 static functype func[32]; 157 static bool init = false; 158 if(!init) 159 { 160 for(int k=0; k<32; ++k) 161 func[k] = 0; 162 163 func[NOTR | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0, true, ColMajor,false,ColMajor,false,ColMajor>::run); 164 func[TR | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0, true, RowMajor,false,ColMajor,false,ColMajor>::run); 165 func[ADJ | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0, true, RowMajor,Conj, ColMajor,false,ColMajor>::run); 166 167 func[NOTR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0, false,ColMajor,false,ColMajor,false,ColMajor>::run); 168 func[TR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0, false,ColMajor,false,RowMajor,false,ColMajor>::run); 169 func[ADJ | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0, false,ColMajor,false,RowMajor,Conj, ColMajor>::run); 170 171 func[NOTR | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0, true, ColMajor,false,ColMajor,false,ColMajor>::run); 172 func[TR | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0, true, RowMajor,false,ColMajor,false,ColMajor>::run); 173 func[ADJ | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0, true, RowMajor,Conj, ColMajor,false,ColMajor>::run); 174 175 func[NOTR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0, false,ColMajor,false,ColMajor,false,ColMajor>::run); 176 func[TR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0, false,ColMajor,false,RowMajor,false,ColMajor>::run); 177 func[ADJ | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0, false,ColMajor,false,RowMajor,Conj, ColMajor>::run); 178 179 func[NOTR | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,true, ColMajor,false,ColMajor,false,ColMajor>::run); 180 func[TR | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,true, RowMajor,false,ColMajor,false,ColMajor>::run); 181 func[ADJ | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,true, RowMajor,Conj, ColMajor,false,ColMajor>::run); 182 183 func[NOTR | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,false,ColMajor,false,ColMajor,false,ColMajor>::run); 184 func[TR | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,false,ColMajor,false,RowMajor,false,ColMajor>::run); 185 func[ADJ | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,false,ColMajor,false,RowMajor,Conj, ColMajor>::run); 186 187 func[NOTR | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,true, ColMajor,false,ColMajor,false,ColMajor>::run); 188 func[TR | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,true, RowMajor,false,ColMajor,false,ColMajor>::run); 189 func[ADJ | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,true, RowMajor,Conj, ColMajor,false,ColMajor>::run); 190 191 func[NOTR | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,false,ColMajor,false,ColMajor,false,ColMajor>::run); 192 func[TR | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,false,ColMajor,false,RowMajor,false,ColMajor>::run); 193 func[ADJ | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,false,ColMajor,false,RowMajor,Conj, ColMajor>::run); 194 195 init = true; 196 } 197 198 Scalar* a = reinterpret_cast<Scalar*>(pa); 199 Scalar* b = reinterpret_cast<Scalar*>(pb); 200 Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 201 202 int info = 0; 203 if(SIDE(*side)==INVALID) info = 1; 204 else if(UPLO(*uplo)==INVALID) info = 2; 205 else if(OP(*opa)==INVALID) info = 3; 206 else if(DIAG(*diag)==INVALID) info = 4; 207 else if(*m<0) info = 5; 208 else if(*n<0) info = 6; 209 else if(*lda<std::max(1,(SIDE(*side)==LEFT)?*m:*n)) info = 9; 210 else if(*ldb<std::max(1,*m)) info = 11; 211 if(info) 212 return xerbla_(SCALAR_SUFFIX_UP"TRMM ",&info,6); 213 214 int code = OP(*opa) | (SIDE(*side) << 2) | (UPLO(*uplo) << 3) | (DIAG(*diag) << 4); 215 216 if(*m==0 || *n==0) 217 return 1; 218 219 // FIXME find a way to avoid this copy 220 Matrix<Scalar,Dynamic,Dynamic,ColMajor> tmp = matrix(b,*m,*n,*ldb); 221 matrix(b,*m,*n,*ldb).setZero(); 222 223 if(SIDE(*side)==LEFT) 224 { 225 internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic,4> blocking(*m,*n,*m); 226 func[code](*m, *n, *m, a, *lda, tmp.data(), tmp.outerStride(), b, *ldb, alpha, blocking); 227 } 228 else 229 { 230 internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic,4> blocking(*m,*n,*n); 231 func[code](*m, *n, *n, tmp.data(), tmp.outerStride(), a, *lda, b, *ldb, alpha, blocking); 232 } 233 return 1; 234} 235 236// c = alpha*a*b + beta*c for side = 'L'or'l' 237// c = alpha*b*a + beta*c for side = 'R'or'r 238int EIGEN_BLAS_FUNC(symm)(char *side, char *uplo, int *m, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc) 239{ 240// std::cerr << "in symm " << *side << " " << *uplo << " " << *m << "x" << *n << " lda:" << *lda << " ldb:" << *ldb << " ldc:" << *ldc << " alpha:" << *palpha << " beta:" << *pbeta << "\n"; 241 Scalar* a = reinterpret_cast<Scalar*>(pa); 242 Scalar* b = reinterpret_cast<Scalar*>(pb); 243 Scalar* c = reinterpret_cast<Scalar*>(pc); 244 Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 245 Scalar beta = *reinterpret_cast<Scalar*>(pbeta); 246 247 int info = 0; 248 if(SIDE(*side)==INVALID) info = 1; 249 else if(UPLO(*uplo)==INVALID) info = 2; 250 else if(*m<0) info = 3; 251 else if(*n<0) info = 4; 252 else if(*lda<std::max(1,(SIDE(*side)==LEFT)?*m:*n)) info = 7; 253 else if(*ldb<std::max(1,*m)) info = 9; 254 else if(*ldc<std::max(1,*m)) info = 12; 255 if(info) 256 return xerbla_(SCALAR_SUFFIX_UP"SYMM ",&info,6); 257 258 if(beta!=Scalar(1)) 259 { 260 if(beta==Scalar(0)) matrix(c, *m, *n, *ldc).setZero(); 261 else matrix(c, *m, *n, *ldc) *= beta; 262 } 263 264 if(*m==0 || *n==0) 265 { 266 return 1; 267 } 268 269 #if ISCOMPLEX 270 // FIXME add support for symmetric complex matrix 271 int size = (SIDE(*side)==LEFT) ? (*m) : (*n); 272 Matrix<Scalar,Dynamic,Dynamic,ColMajor> matA(size,size); 273 if(UPLO(*uplo)==UP) 274 { 275 matA.triangularView<Upper>() = matrix(a,size,size,*lda); 276 matA.triangularView<Lower>() = matrix(a,size,size,*lda).transpose(); 277 } 278 else if(UPLO(*uplo)==LO) 279 { 280 matA.triangularView<Lower>() = matrix(a,size,size,*lda); 281 matA.triangularView<Upper>() = matrix(a,size,size,*lda).transpose(); 282 } 283 if(SIDE(*side)==LEFT) 284 matrix(c, *m, *n, *ldc) += alpha * matA * matrix(b, *m, *n, *ldb); 285 else if(SIDE(*side)==RIGHT) 286 matrix(c, *m, *n, *ldc) += alpha * matrix(b, *m, *n, *ldb) * matA; 287 #else 288 if(SIDE(*side)==LEFT) 289 if(UPLO(*uplo)==UP) internal::product_selfadjoint_matrix<Scalar, DenseIndex, RowMajor,true,false, ColMajor,false,false, ColMajor>::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha); 290 else if(UPLO(*uplo)==LO) internal::product_selfadjoint_matrix<Scalar, DenseIndex, ColMajor,true,false, ColMajor,false,false, ColMajor>::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha); 291 else return 0; 292 else if(SIDE(*side)==RIGHT) 293 if(UPLO(*uplo)==UP) internal::product_selfadjoint_matrix<Scalar, DenseIndex, ColMajor,false,false, RowMajor,true,false, ColMajor>::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha); 294 else if(UPLO(*uplo)==LO) internal::product_selfadjoint_matrix<Scalar, DenseIndex, ColMajor,false,false, ColMajor,true,false, ColMajor>::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha); 295 else return 0; 296 else 297 return 0; 298 #endif 299 300 return 0; 301} 302 303// c = alpha*a*a' + beta*c for op = 'N'or'n' 304// c = alpha*a'*a + beta*c for op = 'T'or't','C'or'c' 305int EIGEN_BLAS_FUNC(syrk)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pbeta, RealScalar *pc, int *ldc) 306{ 307// std::cerr << "in syrk " << *uplo << " " << *op << " " << *n << " " << *k << " " << *palpha << " " << *lda << " " << *pbeta << " " << *ldc << "\n"; 308 #if !ISCOMPLEX 309 typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, const Scalar&); 310 static functype func[8]; 311 312 static bool init = false; 313 if(!init) 314 { 315 for(int k=0; k<8; ++k) 316 func[k] = 0; 317 318 func[NOTR | (UP << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,ColMajor,Conj, Upper>::run); 319 func[TR | (UP << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,false,Scalar,ColMajor,ColMajor,Conj, Upper>::run); 320 func[ADJ | (UP << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,ColMajor,false,Upper>::run); 321 322 func[NOTR | (LO << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,ColMajor,Conj, Lower>::run); 323 func[TR | (LO << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,false,Scalar,ColMajor,ColMajor,Conj, Lower>::run); 324 func[ADJ | (LO << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,ColMajor,false,Lower>::run); 325 326 init = true; 327 } 328 #endif 329 330 Scalar* a = reinterpret_cast<Scalar*>(pa); 331 Scalar* c = reinterpret_cast<Scalar*>(pc); 332 Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 333 Scalar beta = *reinterpret_cast<Scalar*>(pbeta); 334 335 int info = 0; 336 if(UPLO(*uplo)==INVALID) info = 1; 337 else if(OP(*op)==INVALID) info = 2; 338 else if(*n<0) info = 3; 339 else if(*k<0) info = 4; 340 else if(*lda<std::max(1,(OP(*op)==NOTR)?*n:*k)) info = 7; 341 else if(*ldc<std::max(1,*n)) info = 10; 342 if(info) 343 return xerbla_(SCALAR_SUFFIX_UP"SYRK ",&info,6); 344 345 if(beta!=Scalar(1)) 346 { 347 if(UPLO(*uplo)==UP) 348 if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero(); 349 else matrix(c, *n, *n, *ldc).triangularView<Upper>() *= beta; 350 else 351 if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero(); 352 else matrix(c, *n, *n, *ldc).triangularView<Lower>() *= beta; 353 } 354 355 #if ISCOMPLEX 356 // FIXME add support for symmetric complex matrix 357 if(UPLO(*uplo)==UP) 358 { 359 if(OP(*op)==NOTR) 360 matrix(c, *n, *n, *ldc).triangularView<Upper>() += alpha * matrix(a,*n,*k,*lda) * matrix(a,*n,*k,*lda).transpose(); 361 else 362 matrix(c, *n, *n, *ldc).triangularView<Upper>() += alpha * matrix(a,*k,*n,*lda).transpose() * matrix(a,*k,*n,*lda); 363 } 364 else 365 { 366 if(OP(*op)==NOTR) 367 matrix(c, *n, *n, *ldc).triangularView<Lower>() += alpha * matrix(a,*n,*k,*lda) * matrix(a,*n,*k,*lda).transpose(); 368 else 369 matrix(c, *n, *n, *ldc).triangularView<Lower>() += alpha * matrix(a,*k,*n,*lda).transpose() * matrix(a,*k,*n,*lda); 370 } 371 #else 372 int code = OP(*op) | (UPLO(*uplo) << 2); 373 func[code](*n, *k, a, *lda, a, *lda, c, *ldc, alpha); 374 #endif 375 376 return 0; 377} 378 379// c = alpha*a*b' + alpha*b*a' + beta*c for op = 'N'or'n' 380// c = alpha*a'*b + alpha*b'*a + beta*c for op = 'T'or't' 381int EIGEN_BLAS_FUNC(syr2k)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc) 382{ 383 Scalar* a = reinterpret_cast<Scalar*>(pa); 384 Scalar* b = reinterpret_cast<Scalar*>(pb); 385 Scalar* c = reinterpret_cast<Scalar*>(pc); 386 Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 387 Scalar beta = *reinterpret_cast<Scalar*>(pbeta); 388 389 int info = 0; 390 if(UPLO(*uplo)==INVALID) info = 1; 391 else if(OP(*op)==INVALID) info = 2; 392 else if(*n<0) info = 3; 393 else if(*k<0) info = 4; 394 else if(*lda<std::max(1,(OP(*op)==NOTR)?*n:*k)) info = 7; 395 else if(*ldb<std::max(1,(OP(*op)==NOTR)?*n:*k)) info = 9; 396 else if(*ldc<std::max(1,*n)) info = 12; 397 if(info) 398 return xerbla_(SCALAR_SUFFIX_UP"SYR2K",&info,6); 399 400 if(beta!=Scalar(1)) 401 { 402 if(UPLO(*uplo)==UP) 403 if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero(); 404 else matrix(c, *n, *n, *ldc).triangularView<Upper>() *= beta; 405 else 406 if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero(); 407 else matrix(c, *n, *n, *ldc).triangularView<Lower>() *= beta; 408 } 409 410 if(*k==0) 411 return 1; 412 413 if(OP(*op)==NOTR) 414 { 415 if(UPLO(*uplo)==UP) 416 { 417 matrix(c, *n, *n, *ldc).triangularView<Upper>() 418 += alpha *matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).transpose() 419 + alpha*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).transpose(); 420 } 421 else if(UPLO(*uplo)==LO) 422 matrix(c, *n, *n, *ldc).triangularView<Lower>() 423 += alpha*matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).transpose() 424 + alpha*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).transpose(); 425 } 426 else if(OP(*op)==TR || OP(*op)==ADJ) 427 { 428 if(UPLO(*uplo)==UP) 429 matrix(c, *n, *n, *ldc).triangularView<Upper>() 430 += alpha*matrix(a, *k, *n, *lda).transpose()*matrix(b, *k, *n, *ldb) 431 + alpha*matrix(b, *k, *n, *ldb).transpose()*matrix(a, *k, *n, *lda); 432 else if(UPLO(*uplo)==LO) 433 matrix(c, *n, *n, *ldc).triangularView<Lower>() 434 += alpha*matrix(a, *k, *n, *lda).transpose()*matrix(b, *k, *n, *ldb) 435 + alpha*matrix(b, *k, *n, *ldb).transpose()*matrix(a, *k, *n, *lda); 436 } 437 438 return 0; 439} 440 441 442#if ISCOMPLEX 443 444// c = alpha*a*b + beta*c for side = 'L'or'l' 445// c = alpha*b*a + beta*c for side = 'R'or'r 446int EIGEN_BLAS_FUNC(hemm)(char *side, char *uplo, int *m, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc) 447{ 448 Scalar* a = reinterpret_cast<Scalar*>(pa); 449 Scalar* b = reinterpret_cast<Scalar*>(pb); 450 Scalar* c = reinterpret_cast<Scalar*>(pc); 451 Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 452 Scalar beta = *reinterpret_cast<Scalar*>(pbeta); 453 454// std::cerr << "in hemm " << *side << " " << *uplo << " " << *m << " " << *n << " " << alpha << " " << *lda << " " << beta << " " << *ldc << "\n"; 455 456 int info = 0; 457 if(SIDE(*side)==INVALID) info = 1; 458 else if(UPLO(*uplo)==INVALID) info = 2; 459 else if(*m<0) info = 3; 460 else if(*n<0) info = 4; 461 else if(*lda<std::max(1,(SIDE(*side)==LEFT)?*m:*n)) info = 7; 462 else if(*ldb<std::max(1,*m)) info = 9; 463 else if(*ldc<std::max(1,*m)) info = 12; 464 if(info) 465 return xerbla_(SCALAR_SUFFIX_UP"HEMM ",&info,6); 466 467 if(beta==Scalar(0)) matrix(c, *m, *n, *ldc).setZero(); 468 else if(beta!=Scalar(1)) matrix(c, *m, *n, *ldc) *= beta; 469 470 if(*m==0 || *n==0) 471 { 472 return 1; 473 } 474 475 if(SIDE(*side)==LEFT) 476 { 477 if(UPLO(*uplo)==UP) internal::product_selfadjoint_matrix<Scalar,DenseIndex,RowMajor,true,Conj, ColMajor,false,false, ColMajor> 478 ::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha); 479 else if(UPLO(*uplo)==LO) internal::product_selfadjoint_matrix<Scalar,DenseIndex,ColMajor,true,false, ColMajor,false,false, ColMajor> 480 ::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha); 481 else return 0; 482 } 483 else if(SIDE(*side)==RIGHT) 484 { 485 if(UPLO(*uplo)==UP) matrix(c,*m,*n,*ldc) += alpha * matrix(b,*m,*n,*ldb) * matrix(a,*n,*n,*lda).selfadjointView<Upper>();/*internal::product_selfadjoint_matrix<Scalar,DenseIndex,ColMajor,false,false, RowMajor,true,Conj, ColMajor> 486 ::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha);*/ 487 else if(UPLO(*uplo)==LO) internal::product_selfadjoint_matrix<Scalar,DenseIndex,ColMajor,false,false, ColMajor,true,false, ColMajor> 488 ::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha); 489 else return 0; 490 } 491 else 492 { 493 return 0; 494 } 495 496 return 0; 497} 498 499// c = alpha*a*conj(a') + beta*c for op = 'N'or'n' 500// c = alpha*conj(a')*a + beta*c for op = 'C'or'c' 501int EIGEN_BLAS_FUNC(herk)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pbeta, RealScalar *pc, int *ldc) 502{ 503 typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, const Scalar&); 504 static functype func[8]; 505 506 static bool init = false; 507 if(!init) 508 { 509 for(int k=0; k<8; ++k) 510 func[k] = 0; 511 512 func[NOTR | (UP << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,Conj, ColMajor,Upper>::run); 513 func[ADJ | (UP << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,false,ColMajor,Upper>::run); 514 515 func[NOTR | (LO << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,Conj, ColMajor,Lower>::run); 516 func[ADJ | (LO << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,false,ColMajor,Lower>::run); 517 518 init = true; 519 } 520 521 Scalar* a = reinterpret_cast<Scalar*>(pa); 522 Scalar* c = reinterpret_cast<Scalar*>(pc); 523 RealScalar alpha = *palpha; 524 RealScalar beta = *pbeta; 525 526// std::cerr << "in herk " << *uplo << " " << *op << " " << *n << " " << *k << " " << alpha << " " << *lda << " " << beta << " " << *ldc << "\n"; 527 528 int info = 0; 529 if(UPLO(*uplo)==INVALID) info = 1; 530 else if((OP(*op)==INVALID) || (OP(*op)==TR)) info = 2; 531 else if(*n<0) info = 3; 532 else if(*k<0) info = 4; 533 else if(*lda<std::max(1,(OP(*op)==NOTR)?*n:*k)) info = 7; 534 else if(*ldc<std::max(1,*n)) info = 10; 535 if(info) 536 return xerbla_(SCALAR_SUFFIX_UP"HERK ",&info,6); 537 538 int code = OP(*op) | (UPLO(*uplo) << 2); 539 540 if(beta!=RealScalar(1)) 541 { 542 if(UPLO(*uplo)==UP) 543 if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero(); 544 else matrix(c, *n, *n, *ldc).triangularView<StrictlyUpper>() *= beta; 545 else 546 if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero(); 547 else matrix(c, *n, *n, *ldc).triangularView<StrictlyLower>() *= beta; 548 549 if(beta!=Scalar(0)) 550 { 551 matrix(c, *n, *n, *ldc).diagonal().real() *= beta; 552 matrix(c, *n, *n, *ldc).diagonal().imag().setZero(); 553 } 554 } 555 556 if(*k>0 && alpha!=RealScalar(0)) 557 { 558 func[code](*n, *k, a, *lda, a, *lda, c, *ldc, alpha); 559 matrix(c, *n, *n, *ldc).diagonal().imag().setZero(); 560 } 561 return 0; 562} 563 564// c = alpha*a*conj(b') + conj(alpha)*b*conj(a') + beta*c, for op = 'N'or'n' 565// c = alpha*conj(a')*b + conj(alpha)*conj(b')*a + beta*c, for op = 'C'or'c' 566int EIGEN_BLAS_FUNC(her2k)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc) 567{ 568 Scalar* a = reinterpret_cast<Scalar*>(pa); 569 Scalar* b = reinterpret_cast<Scalar*>(pb); 570 Scalar* c = reinterpret_cast<Scalar*>(pc); 571 Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 572 RealScalar beta = *pbeta; 573 574 int info = 0; 575 if(UPLO(*uplo)==INVALID) info = 1; 576 else if((OP(*op)==INVALID) || (OP(*op)==TR)) info = 2; 577 else if(*n<0) info = 3; 578 else if(*k<0) info = 4; 579 else if(*lda<std::max(1,(OP(*op)==NOTR)?*n:*k)) info = 7; 580 else if(*lda<std::max(1,(OP(*op)==NOTR)?*n:*k)) info = 9; 581 else if(*ldc<std::max(1,*n)) info = 12; 582 if(info) 583 return xerbla_(SCALAR_SUFFIX_UP"HER2K",&info,6); 584 585 if(beta!=RealScalar(1)) 586 { 587 if(UPLO(*uplo)==UP) 588 if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero(); 589 else matrix(c, *n, *n, *ldc).triangularView<StrictlyUpper>() *= beta; 590 else 591 if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero(); 592 else matrix(c, *n, *n, *ldc).triangularView<StrictlyLower>() *= beta; 593 594 if(beta!=Scalar(0)) 595 { 596 matrix(c, *n, *n, *ldc).diagonal().real() *= beta; 597 matrix(c, *n, *n, *ldc).diagonal().imag().setZero(); 598 } 599 } 600 else if(*k>0 && alpha!=Scalar(0)) 601 matrix(c, *n, *n, *ldc).diagonal().imag().setZero(); 602 603 if(*k==0) 604 return 1; 605 606 if(OP(*op)==NOTR) 607 { 608 if(UPLO(*uplo)==UP) 609 { 610 matrix(c, *n, *n, *ldc).triangularView<Upper>() 611 += alpha *matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).adjoint() 612 + numext::conj(alpha)*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).adjoint(); 613 } 614 else if(UPLO(*uplo)==LO) 615 matrix(c, *n, *n, *ldc).triangularView<Lower>() 616 += alpha*matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).adjoint() 617 + numext::conj(alpha)*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).adjoint(); 618 } 619 else if(OP(*op)==ADJ) 620 { 621 if(UPLO(*uplo)==UP) 622 matrix(c, *n, *n, *ldc).triangularView<Upper>() 623 += alpha*matrix(a, *k, *n, *lda).adjoint()*matrix(b, *k, *n, *ldb) 624 + numext::conj(alpha)*matrix(b, *k, *n, *ldb).adjoint()*matrix(a, *k, *n, *lda); 625 else if(UPLO(*uplo)==LO) 626 matrix(c, *n, *n, *ldc).triangularView<Lower>() 627 += alpha*matrix(a, *k, *n, *lda).adjoint()*matrix(b, *k, *n, *ldb) 628 + numext::conj(alpha)*matrix(b, *k, *n, *ldb).adjoint()*matrix(a, *k, *n, *lda); 629 } 630 631 return 1; 632} 633 634#endif // ISCOMPLEX 635