• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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 
EIGEN_BLAS_FUNC(gemm)12 int 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 
EIGEN_BLAS_FUNC(trsm)66 int 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'
EIGEN_BLAS_FUNC(trmm)152 int 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, 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
EIGEN_BLAS_FUNC(symm)238 int 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'
EIGEN_BLAS_FUNC(syrk)305 int 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   typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, Scalar);
309   static functype func[8];
310 
311   static bool init = false;
312   if(!init)
313   {
314     for(int k=0; k<8; ++k)
315       func[k] = 0;
316 
317     func[NOTR  | (UP << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,ColMajor,Conj, Upper>::run);
318     func[TR    | (UP << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,false,Scalar,ColMajor,ColMajor,Conj, Upper>::run);
319     func[ADJ   | (UP << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,ColMajor,false,Upper>::run);
320 
321     func[NOTR  | (LO << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,ColMajor,Conj, Lower>::run);
322     func[TR    | (LO << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,false,Scalar,ColMajor,ColMajor,Conj, Lower>::run);
323     func[ADJ   | (LO << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,ColMajor,false,Lower>::run);
324 
325     init = true;
326   }
327 
328   Scalar* a = reinterpret_cast<Scalar*>(pa);
329   Scalar* c = reinterpret_cast<Scalar*>(pc);
330   Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
331   Scalar beta  = *reinterpret_cast<Scalar*>(pbeta);
332 
333   int info = 0;
334   if(UPLO(*uplo)==INVALID)                                            info = 1;
335   else if(OP(*op)==INVALID)                                           info = 2;
336   else if(*n<0)                                                       info = 3;
337   else if(*k<0)                                                       info = 4;
338   else if(*lda<std::max(1,(OP(*op)==NOTR)?*n:*k))                     info = 7;
339   else if(*ldc<std::max(1,*n))                                        info = 10;
340   if(info)
341     return xerbla_(SCALAR_SUFFIX_UP"SYRK ",&info,6);
342 
343   if(beta!=Scalar(1))
344   {
345     if(UPLO(*uplo)==UP)
346       if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero();
347       else                matrix(c, *n, *n, *ldc).triangularView<Upper>() *= beta;
348     else
349       if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero();
350       else                matrix(c, *n, *n, *ldc).triangularView<Lower>() *= beta;
351   }
352 
353   #if ISCOMPLEX
354   // FIXME add support for symmetric complex matrix
355   if(UPLO(*uplo)==UP)
356   {
357     if(OP(*op)==NOTR)
358       matrix(c, *n, *n, *ldc).triangularView<Upper>() += alpha * matrix(a,*n,*k,*lda) * matrix(a,*n,*k,*lda).transpose();
359     else
360       matrix(c, *n, *n, *ldc).triangularView<Upper>() += alpha * matrix(a,*k,*n,*lda).transpose() * matrix(a,*k,*n,*lda);
361   }
362   else
363   {
364     if(OP(*op)==NOTR)
365       matrix(c, *n, *n, *ldc).triangularView<Lower>() += alpha * matrix(a,*n,*k,*lda) * matrix(a,*n,*k,*lda).transpose();
366     else
367       matrix(c, *n, *n, *ldc).triangularView<Lower>() += alpha * matrix(a,*k,*n,*lda).transpose() * matrix(a,*k,*n,*lda);
368   }
369   #else
370   int code = OP(*op) | (UPLO(*uplo) << 2);
371   func[code](*n, *k, a, *lda, a, *lda, c, *ldc, alpha);
372   #endif
373 
374   return 0;
375 }
376 
377 // c = alpha*a*b' + alpha*b*a' + beta*c  for op = 'N'or'n'
378 // c = alpha*a'*b + alpha*b'*a + beta*c  for op = 'T'or't'
EIGEN_BLAS_FUNC(syr2k)379 int 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)
380 {
381   Scalar* a = reinterpret_cast<Scalar*>(pa);
382   Scalar* b = reinterpret_cast<Scalar*>(pb);
383   Scalar* c = reinterpret_cast<Scalar*>(pc);
384   Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
385   Scalar beta  = *reinterpret_cast<Scalar*>(pbeta);
386 
387   int info = 0;
388   if(UPLO(*uplo)==INVALID)                                            info = 1;
389   else if(OP(*op)==INVALID)                                           info = 2;
390   else if(*n<0)                                                       info = 3;
391   else if(*k<0)                                                       info = 4;
392   else if(*lda<std::max(1,(OP(*op)==NOTR)?*n:*k))                     info = 7;
393   else if(*ldb<std::max(1,(OP(*op)==NOTR)?*n:*k))                     info = 9;
394   else if(*ldc<std::max(1,*n))                                        info = 12;
395   if(info)
396     return xerbla_(SCALAR_SUFFIX_UP"SYR2K",&info,6);
397 
398   if(beta!=Scalar(1))
399   {
400     if(UPLO(*uplo)==UP)
401       if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero();
402       else                matrix(c, *n, *n, *ldc).triangularView<Upper>() *= beta;
403     else
404       if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero();
405       else                matrix(c, *n, *n, *ldc).triangularView<Lower>() *= beta;
406   }
407 
408   if(*k==0)
409     return 1;
410 
411   if(OP(*op)==NOTR)
412   {
413     if(UPLO(*uplo)==UP)
414     {
415       matrix(c, *n, *n, *ldc).triangularView<Upper>()
416         += alpha *matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).transpose()
417         +  alpha*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).transpose();
418     }
419     else if(UPLO(*uplo)==LO)
420       matrix(c, *n, *n, *ldc).triangularView<Lower>()
421         += alpha*matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).transpose()
422         +  alpha*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).transpose();
423   }
424   else if(OP(*op)==TR || OP(*op)==ADJ)
425   {
426     if(UPLO(*uplo)==UP)
427       matrix(c, *n, *n, *ldc).triangularView<Upper>()
428         += alpha*matrix(a, *k, *n, *lda).transpose()*matrix(b, *k, *n, *ldb)
429         +  alpha*matrix(b, *k, *n, *ldb).transpose()*matrix(a, *k, *n, *lda);
430     else if(UPLO(*uplo)==LO)
431       matrix(c, *n, *n, *ldc).triangularView<Lower>()
432         += alpha*matrix(a, *k, *n, *lda).transpose()*matrix(b, *k, *n, *ldb)
433         +  alpha*matrix(b, *k, *n, *ldb).transpose()*matrix(a, *k, *n, *lda);
434   }
435 
436   return 0;
437 }
438 
439 
440 #if ISCOMPLEX
441 
442 // c = alpha*a*b + beta*c  for side = 'L'or'l'
443 // c = alpha*b*a + beta*c  for side = 'R'or'r
EIGEN_BLAS_FUNC(hemm)444 int 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)
445 {
446   Scalar* a = reinterpret_cast<Scalar*>(pa);
447   Scalar* b = reinterpret_cast<Scalar*>(pb);
448   Scalar* c = reinterpret_cast<Scalar*>(pc);
449   Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
450   Scalar beta  = *reinterpret_cast<Scalar*>(pbeta);
451 
452 //   std::cerr << "in hemm " << *side << " " << *uplo << " " << *m << " " << *n << " " << alpha << " " << *lda << " " << beta << " " << *ldc << "\n";
453 
454   int info = 0;
455   if(SIDE(*side)==INVALID)                                            info = 1;
456   else if(UPLO(*uplo)==INVALID)                                       info = 2;
457   else if(*m<0)                                                       info = 3;
458   else if(*n<0)                                                       info = 4;
459   else if(*lda<std::max(1,(SIDE(*side)==LEFT)?*m:*n))                 info = 7;
460   else if(*ldb<std::max(1,*m))                                        info = 9;
461   else if(*ldc<std::max(1,*m))                                        info = 12;
462   if(info)
463     return xerbla_(SCALAR_SUFFIX_UP"HEMM ",&info,6);
464 
465   if(beta==Scalar(0))       matrix(c, *m, *n, *ldc).setZero();
466   else if(beta!=Scalar(1))  matrix(c, *m, *n, *ldc) *= beta;
467 
468   if(*m==0 || *n==0)
469   {
470     return 1;
471   }
472 
473   if(SIDE(*side)==LEFT)
474   {
475     if(UPLO(*uplo)==UP)       internal::product_selfadjoint_matrix<Scalar,DenseIndex,RowMajor,true,Conj,  ColMajor,false,false, ColMajor>
476                                 ::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha);
477     else if(UPLO(*uplo)==LO)  internal::product_selfadjoint_matrix<Scalar,DenseIndex,ColMajor,true,false, ColMajor,false,false, ColMajor>
478                                 ::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha);
479     else                      return 0;
480   }
481   else if(SIDE(*side)==RIGHT)
482   {
483     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>
484                                 ::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha);*/
485     else if(UPLO(*uplo)==LO)  internal::product_selfadjoint_matrix<Scalar,DenseIndex,ColMajor,false,false, ColMajor,true,false, ColMajor>
486                                 ::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha);
487     else                      return 0;
488   }
489   else
490   {
491     return 0;
492   }
493 
494   return 0;
495 }
496 
497 // c = alpha*a*conj(a') + beta*c  for op = 'N'or'n'
498 // c = alpha*conj(a')*a + beta*c  for op  = 'C'or'c'
EIGEN_BLAS_FUNC(herk)499 int EIGEN_BLAS_FUNC(herk)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pbeta, RealScalar *pc, int *ldc)
500 {
501   typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, Scalar);
502   static functype func[8];
503 
504   static bool init = false;
505   if(!init)
506   {
507     for(int k=0; k<8; ++k)
508       func[k] = 0;
509 
510     func[NOTR  | (UP << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,Conj, ColMajor,Upper>::run);
511     func[ADJ   | (UP << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,false,ColMajor,Upper>::run);
512 
513     func[NOTR  | (LO << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,Conj, ColMajor,Lower>::run);
514     func[ADJ   | (LO << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,false,ColMajor,Lower>::run);
515 
516     init = true;
517   }
518 
519   Scalar* a = reinterpret_cast<Scalar*>(pa);
520   Scalar* c = reinterpret_cast<Scalar*>(pc);
521   RealScalar alpha = *palpha;
522   RealScalar beta  = *pbeta;
523 
524 //   std::cerr << "in herk " << *uplo << " " << *op << " " << *n << " " << *k << " " << alpha << " " << *lda << " " << beta << " " << *ldc << "\n";
525 
526   int info = 0;
527   if(UPLO(*uplo)==INVALID)                                            info = 1;
528   else if((OP(*op)==INVALID) || (OP(*op)==TR))                        info = 2;
529   else if(*n<0)                                                       info = 3;
530   else if(*k<0)                                                       info = 4;
531   else if(*lda<std::max(1,(OP(*op)==NOTR)?*n:*k))                     info = 7;
532   else if(*ldc<std::max(1,*n))                                        info = 10;
533   if(info)
534     return xerbla_(SCALAR_SUFFIX_UP"HERK ",&info,6);
535 
536   int code = OP(*op) | (UPLO(*uplo) << 2);
537 
538   if(beta!=RealScalar(1))
539   {
540     if(UPLO(*uplo)==UP)
541       if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero();
542       else                matrix(c, *n, *n, *ldc).triangularView<StrictlyUpper>() *= beta;
543     else
544       if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero();
545       else                matrix(c, *n, *n, *ldc).triangularView<StrictlyLower>() *= beta;
546 
547     if(beta!=Scalar(0))
548     {
549       matrix(c, *n, *n, *ldc).diagonal().real() *= beta;
550       matrix(c, *n, *n, *ldc).diagonal().imag().setZero();
551     }
552   }
553 
554   if(*k>0 && alpha!=RealScalar(0))
555   {
556     func[code](*n, *k, a, *lda, a, *lda, c, *ldc, alpha);
557     matrix(c, *n, *n, *ldc).diagonal().imag().setZero();
558   }
559   return 0;
560 }
561 
562 // c = alpha*a*conj(b') + conj(alpha)*b*conj(a') + beta*c,  for op = 'N'or'n'
563 // c = alpha*conj(a')*b + conj(alpha)*conj(b')*a + beta*c,  for op = 'C'or'c'
EIGEN_BLAS_FUNC(her2k)564 int 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)
565 {
566   Scalar* a = reinterpret_cast<Scalar*>(pa);
567   Scalar* b = reinterpret_cast<Scalar*>(pb);
568   Scalar* c = reinterpret_cast<Scalar*>(pc);
569   Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
570   RealScalar beta  = *pbeta;
571 
572   int info = 0;
573   if(UPLO(*uplo)==INVALID)                                            info = 1;
574   else if((OP(*op)==INVALID) || (OP(*op)==TR))                        info = 2;
575   else if(*n<0)                                                       info = 3;
576   else if(*k<0)                                                       info = 4;
577   else if(*lda<std::max(1,(OP(*op)==NOTR)?*n:*k))                     info = 7;
578   else if(*lda<std::max(1,(OP(*op)==NOTR)?*n:*k))                     info = 9;
579   else if(*ldc<std::max(1,*n))                                        info = 12;
580   if(info)
581     return xerbla_(SCALAR_SUFFIX_UP"HER2K",&info,6);
582 
583   if(beta!=RealScalar(1))
584   {
585     if(UPLO(*uplo)==UP)
586       if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero();
587       else                matrix(c, *n, *n, *ldc).triangularView<StrictlyUpper>() *= beta;
588     else
589       if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero();
590       else                matrix(c, *n, *n, *ldc).triangularView<StrictlyLower>() *= beta;
591 
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 
601   if(*k==0)
602     return 1;
603 
604   if(OP(*op)==NOTR)
605   {
606     if(UPLO(*uplo)==UP)
607     {
608       matrix(c, *n, *n, *ldc).triangularView<Upper>()
609         +=         alpha *matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).adjoint()
610         +  internal::conj(alpha)*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).adjoint();
611     }
612     else if(UPLO(*uplo)==LO)
613       matrix(c, *n, *n, *ldc).triangularView<Lower>()
614         += alpha*matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).adjoint()
615         +  internal::conj(alpha)*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).adjoint();
616   }
617   else if(OP(*op)==ADJ)
618   {
619     if(UPLO(*uplo)==UP)
620       matrix(c, *n, *n, *ldc).triangularView<Upper>()
621         += alpha*matrix(a, *k, *n, *lda).adjoint()*matrix(b, *k, *n, *ldb)
622         +  internal::conj(alpha)*matrix(b, *k, *n, *ldb).adjoint()*matrix(a, *k, *n, *lda);
623     else if(UPLO(*uplo)==LO)
624       matrix(c, *n, *n, *ldc).triangularView<Lower>()
625         += alpha*matrix(a, *k, *n, *lda).adjoint()*matrix(b, *k, *n, *ldb)
626         +  internal::conj(alpha)*matrix(b, *k, *n, *ldb).adjoint()*matrix(a, *k, *n, *lda);
627   }
628 
629   return 1;
630 }
631 
632 #endif // ISCOMPLEX
633