level3_impl.h 37 KB

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