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(gemv)12 int EIGEN_BLAS_FUNC(gemv)(char *opa, int *m, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *incb, RealScalar *pbeta, RealScalar *pc, int *incc)
13 {
14   typedef void (*functype)(int, int, const Scalar *, int, const Scalar *, int , Scalar *, int, Scalar);
15   static functype func[4];
16 
17   static bool init = false;
18   if(!init)
19   {
20     for(int k=0; k<4; ++k)
21       func[k] = 0;
22 
23     func[NOTR] = (internal::general_matrix_vector_product<int,Scalar,ColMajor,false,Scalar,false>::run);
24     func[TR  ] = (internal::general_matrix_vector_product<int,Scalar,RowMajor,false,Scalar,false>::run);
25     func[ADJ ] = (internal::general_matrix_vector_product<int,Scalar,RowMajor,Conj, Scalar,false>::run);
26 
27     init = true;
28   }
29 
30   Scalar* a = reinterpret_cast<Scalar*>(pa);
31   Scalar* b = reinterpret_cast<Scalar*>(pb);
32   Scalar* c = reinterpret_cast<Scalar*>(pc);
33   Scalar alpha  = *reinterpret_cast<Scalar*>(palpha);
34   Scalar beta   = *reinterpret_cast<Scalar*>(pbeta);
35 
36   // check arguments
37   int info = 0;
38   if(OP(*opa)==INVALID)           info = 1;
39   else if(*m<0)                   info = 2;
40   else if(*n<0)                   info = 3;
41   else if(*lda<std::max(1,*m))    info = 6;
42   else if(*incb==0)               info = 8;
43   else if(*incc==0)               info = 11;
44   if(info)
45     return xerbla_(SCALAR_SUFFIX_UP"GEMV ",&info,6);
46 
47   if(*m==0 || *n==0 || (alpha==Scalar(0) && beta==Scalar(1)))
48     return 0;
49 
50   int actual_m = *m;
51   int actual_n = *n;
52   int code = OP(*opa);
53   if(code!=NOTR)
54     std::swap(actual_m,actual_n);
55 
56   Scalar* actual_b = get_compact_vector(b,actual_n,*incb);
57   Scalar* actual_c = get_compact_vector(c,actual_m,*incc);
58 
59   if(beta!=Scalar(1))
60   {
61     if(beta==Scalar(0)) vector(actual_c, actual_m).setZero();
62     else                vector(actual_c, actual_m) *= beta;
63   }
64 
65   if(code>=4 || func[code]==0)
66     return 0;
67 
68   func[code](actual_m, actual_n, a, *lda, actual_b, 1, actual_c, 1, alpha);
69 
70   if(actual_b!=b) delete[] actual_b;
71   if(actual_c!=c) delete[] copy_back(actual_c,c,actual_m,*incc);
72 
73   return 1;
74 }
75 
EIGEN_BLAS_FUNC(trsv)76 int EIGEN_BLAS_FUNC(trsv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pa, int *lda, RealScalar *pb, int *incb)
77 {
78   typedef void (*functype)(int, const Scalar *, int, Scalar *);
79   static functype func[16];
80 
81   static bool init = false;
82   if(!init)
83   {
84     for(int k=0; k<16; ++k)
85       func[k] = 0;
86 
87     func[NOTR  | (UP << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0,       false,ColMajor>::run);
88     func[TR    | (UP << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0,       false,RowMajor>::run);
89     func[ADJ   | (UP << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0,       Conj, RowMajor>::run);
90 
91     func[NOTR  | (LO << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0,       false,ColMajor>::run);
92     func[TR    | (LO << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0,       false,RowMajor>::run);
93     func[ADJ   | (LO << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0,       Conj, RowMajor>::run);
94 
95     func[NOTR  | (UP << 2) | (UNIT  << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,ColMajor>::run);
96     func[TR    | (UP << 2) | (UNIT  << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,RowMajor>::run);
97     func[ADJ   | (UP << 2) | (UNIT  << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,Conj, RowMajor>::run);
98 
99     func[NOTR  | (LO << 2) | (UNIT  << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,ColMajor>::run);
100     func[TR    | (LO << 2) | (UNIT  << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,RowMajor>::run);
101     func[ADJ   | (LO << 2) | (UNIT  << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,Conj, RowMajor>::run);
102 
103     init = true;
104   }
105 
106   Scalar* a = reinterpret_cast<Scalar*>(pa);
107   Scalar* b = reinterpret_cast<Scalar*>(pb);
108 
109   int info = 0;
110   if(UPLO(*uplo)==INVALID)                                            info = 1;
111   else if(OP(*opa)==INVALID)                                          info = 2;
112   else if(DIAG(*diag)==INVALID)                                       info = 3;
113   else if(*n<0)                                                       info = 4;
114   else if(*lda<std::max(1,*n))                                        info = 6;
115   else if(*incb==0)                                                   info = 8;
116   if(info)
117     return xerbla_(SCALAR_SUFFIX_UP"TRSV ",&info,6);
118 
119   Scalar* actual_b = get_compact_vector(b,*n,*incb);
120 
121   int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
122   func[code](*n, a, *lda, actual_b);
123 
124   if(actual_b!=b) delete[] copy_back(actual_b,b,*n,*incb);
125 
126   return 0;
127 }
128 
129 
130 
EIGEN_BLAS_FUNC(trmv)131 int EIGEN_BLAS_FUNC(trmv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pa, int *lda, RealScalar *pb, int *incb)
132 {
133   typedef void (*functype)(int, int, const Scalar *, int, const Scalar *, int, Scalar *, int, const Scalar&);
134   static functype func[16];
135 
136   static bool init = false;
137   if(!init)
138   {
139     for(int k=0; k<16; ++k)
140       func[k] = 0;
141 
142     func[NOTR  | (UP << 2) | (NUNIT << 3)] = (internal::triangular_matrix_vector_product<int,Upper|0,       Scalar,false,Scalar,false,ColMajor>::run);
143     func[TR    | (UP << 2) | (NUNIT << 3)] = (internal::triangular_matrix_vector_product<int,Lower|0,       Scalar,false,Scalar,false,RowMajor>::run);
144     func[ADJ   | (UP << 2) | (NUNIT << 3)] = (internal::triangular_matrix_vector_product<int,Lower|0,       Scalar,Conj, Scalar,false,RowMajor>::run);
145 
146     func[NOTR  | (LO << 2) | (NUNIT << 3)] = (internal::triangular_matrix_vector_product<int,Lower|0,       Scalar,false,Scalar,false,ColMajor>::run);
147     func[TR    | (LO << 2) | (NUNIT << 3)] = (internal::triangular_matrix_vector_product<int,Upper|0,       Scalar,false,Scalar,false,RowMajor>::run);
148     func[ADJ   | (LO << 2) | (NUNIT << 3)] = (internal::triangular_matrix_vector_product<int,Upper|0,       Scalar,Conj, Scalar,false,RowMajor>::run);
149 
150     func[NOTR  | (UP << 2) | (UNIT  << 3)] = (internal::triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run);
151     func[TR    | (UP << 2) | (UNIT  << 3)] = (internal::triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run);
152     func[ADJ   | (UP << 2) | (UNIT  << 3)] = (internal::triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run);
153 
154     func[NOTR  | (LO << 2) | (UNIT  << 3)] = (internal::triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run);
155     func[TR    | (LO << 2) | (UNIT  << 3)] = (internal::triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run);
156     func[ADJ   | (LO << 2) | (UNIT  << 3)] = (internal::triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run);
157 
158     init = true;
159   }
160 
161   Scalar* a = reinterpret_cast<Scalar*>(pa);
162   Scalar* b = reinterpret_cast<Scalar*>(pb);
163 
164   int info = 0;
165   if(UPLO(*uplo)==INVALID)                                            info = 1;
166   else if(OP(*opa)==INVALID)                                          info = 2;
167   else if(DIAG(*diag)==INVALID)                                       info = 3;
168   else if(*n<0)                                                       info = 4;
169   else if(*lda<std::max(1,*n))                                        info = 6;
170   else if(*incb==0)                                                   info = 8;
171   if(info)
172     return xerbla_(SCALAR_SUFFIX_UP"TRMV ",&info,6);
173 
174   if(*n==0)
175     return 1;
176 
177   Scalar* actual_b = get_compact_vector(b,*n,*incb);
178   Matrix<Scalar,Dynamic,1> res(*n);
179   res.setZero();
180 
181   int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
182   if(code>=16 || func[code]==0)
183     return 0;
184 
185   func[code](*n, *n, a, *lda, actual_b, 1, res.data(), 1, Scalar(1));
186 
187   copy_back(res.data(),b,*n,*incb);
188   if(actual_b!=b) delete[] actual_b;
189 
190   return 1;
191 }
192 
193 /**  GBMV  performs one of the matrix-vector operations
194   *
195   *     y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,
196   *
197   *  where alpha and beta are scalars, x and y are vectors and A is an
198   *  m by n band matrix, with kl sub-diagonals and ku super-diagonals.
199   */
EIGEN_BLAS_FUNC(gbmv)200 int EIGEN_BLAS_FUNC(gbmv)(char *trans, int *m, int *n, int *kl, int *ku, RealScalar *palpha, RealScalar *pa, int *lda,
201                           RealScalar *px, int *incx, RealScalar *pbeta, RealScalar *py, int *incy)
202 {
203   Scalar* a = reinterpret_cast<Scalar*>(pa);
204   Scalar* x = reinterpret_cast<Scalar*>(px);
205   Scalar* y = reinterpret_cast<Scalar*>(py);
206   Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
207   Scalar beta = *reinterpret_cast<Scalar*>(pbeta);
208   int coeff_rows = *kl+*ku+1;
209 
210   int info = 0;
211        if(OP(*trans)==INVALID)                                        info = 1;
212   else if(*m<0)                                                       info = 2;
213   else if(*n<0)                                                       info = 3;
214   else if(*kl<0)                                                      info = 4;
215   else if(*ku<0)                                                      info = 5;
216   else if(*lda<coeff_rows)                                            info = 8;
217   else if(*incx==0)                                                   info = 10;
218   else if(*incy==0)                                                   info = 13;
219   if(info)
220     return xerbla_(SCALAR_SUFFIX_UP"GBMV ",&info,6);
221 
222   if(*m==0 || *n==0 || (alpha==Scalar(0) && beta==Scalar(1)))
223     return 0;
224 
225   int actual_m = *m;
226   int actual_n = *n;
227   if(OP(*trans)!=NOTR)
228     std::swap(actual_m,actual_n);
229 
230   Scalar* actual_x = get_compact_vector(x,actual_n,*incx);
231   Scalar* actual_y = get_compact_vector(y,actual_m,*incy);
232 
233   if(beta!=Scalar(1))
234   {
235     if(beta==Scalar(0)) vector(actual_y, actual_m).setZero();
236     else                vector(actual_y, actual_m) *= beta;
237   }
238 
239   MatrixType mat_coeffs(a,coeff_rows,*n,*lda);
240 
241   int nb = std::min(*n,(*m)+(*ku));
242   for(int j=0; j<nb; ++j)
243   {
244     int start = std::max(0,j - *ku);
245     int end = std::min((*m)-1,j + *kl);
246     int len = end - start + 1;
247     int offset = (*ku) - j + start;
248     if(OP(*trans)==NOTR)
249       vector(actual_y+start,len) += (alpha*actual_x[j]) * mat_coeffs.col(j).segment(offset,len);
250     else if(OP(*trans)==TR)
251       actual_y[j] += alpha * ( mat_coeffs.col(j).segment(offset,len).transpose() * vector(actual_x+start,len) ).value();
252     else
253       actual_y[j] += alpha * ( mat_coeffs.col(j).segment(offset,len).adjoint()   * vector(actual_x+start,len) ).value();
254   }
255 
256   if(actual_x!=x) delete[] actual_x;
257   if(actual_y!=y) delete[] copy_back(actual_y,y,actual_m,*incy);
258 
259   return 0;
260 }
261 
262 #if 0
263 /**  TBMV  performs one of the matrix-vector operations
264   *
265   *     x := A*x,   or   x := A'*x,
266   *
267   *  where x is an n element vector and  A is an n by n unit, or non-unit,
268   *  upper or lower triangular band matrix, with ( k + 1 ) diagonals.
269   */
270 int EIGEN_BLAS_FUNC(tbmv)(char *uplo, char *opa, char *diag, int *n, int *k, RealScalar *pa, int *lda, RealScalar *px, int *incx)
271 {
272   Scalar* a = reinterpret_cast<Scalar*>(pa);
273   Scalar* x = reinterpret_cast<Scalar*>(px);
274   int coeff_rows = *k + 1;
275 
276   int info = 0;
277        if(UPLO(*uplo)==INVALID)                                       info = 1;
278   else if(OP(*opa)==INVALID)                                          info = 2;
279   else if(DIAG(*diag)==INVALID)                                       info = 3;
280   else if(*n<0)                                                       info = 4;
281   else if(*k<0)                                                       info = 5;
282   else if(*lda<coeff_rows)                                            info = 7;
283   else if(*incx==0)                                                   info = 9;
284   if(info)
285     return xerbla_(SCALAR_SUFFIX_UP"TBMV ",&info,6);
286 
287   if(*n==0)
288     return 0;
289 
290   int actual_n = *n;
291 
292   Scalar* actual_x = get_compact_vector(x,actual_n,*incx);
293 
294   MatrixType mat_coeffs(a,coeff_rows,*n,*lda);
295 
296   int ku = UPLO(*uplo)==UPPER ? *k : 0;
297   int kl = UPLO(*uplo)==LOWER ? *k : 0;
298 
299   for(int j=0; j<*n; ++j)
300   {
301     int start = std::max(0,j - ku);
302     int end = std::min((*m)-1,j + kl);
303     int len = end - start + 1;
304     int offset = (ku) - j + start;
305 
306     if(OP(*trans)==NOTR)
307       vector(actual_y+start,len) += (alpha*actual_x[j]) * mat_coeffs.col(j).segment(offset,len);
308     else if(OP(*trans)==TR)
309       actual_y[j] += alpha * ( mat_coeffs.col(j).segment(offset,len).transpose() * vector(actual_x+start,len) ).value();
310     else
311       actual_y[j] += alpha * ( mat_coeffs.col(j).segment(offset,len).adjoint()   * vector(actual_x+start,len) ).value();
312   }
313 
314   if(actual_x!=x) delete[] actual_x;
315   if(actual_y!=y) delete[] copy_back(actual_y,y,actual_m,*incy);
316 
317   return 0;
318 }
319 #endif
320 
321 /**  DTBSV  solves one of the systems of equations
322   *
323   *     A*x = b,   or   A'*x = b,
324   *
325   *  where b and x are n element vectors and A is an n by n unit, or
326   *  non-unit, upper or lower triangular band matrix, with ( k + 1 )
327   *  diagonals.
328   *
329   *  No test for singularity or near-singularity is included in this
330   *  routine. Such tests must be performed before calling this routine.
331   */
EIGEN_BLAS_FUNC(tbsv)332 int EIGEN_BLAS_FUNC(tbsv)(char *uplo, char *op, char *diag, int *n, int *k, RealScalar *pa, int *lda, RealScalar *px, int *incx)
333 {
334   typedef void (*functype)(int, int, const Scalar *, int, Scalar *);
335   static functype func[16];
336 
337   static bool init = false;
338   if(!init)
339   {
340     for(int k=0; k<16; ++k)
341       func[k] = 0;
342 
343     func[NOTR  | (UP << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector<int,Upper|0,       Scalar,false,Scalar,ColMajor>::run);
344     func[TR    | (UP << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector<int,Lower|0,       Scalar,false,Scalar,RowMajor>::run);
345     func[ADJ   | (UP << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector<int,Lower|0,       Scalar,Conj, Scalar,RowMajor>::run);
346 
347     func[NOTR  | (LO << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector<int,Lower|0,       Scalar,false,Scalar,ColMajor>::run);
348     func[TR    | (LO << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector<int,Upper|0,       Scalar,false,Scalar,RowMajor>::run);
349     func[ADJ   | (LO << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector<int,Upper|0,       Scalar,Conj, Scalar,RowMajor>::run);
350 
351     func[NOTR  | (UP << 2) | (UNIT  << 3)] = (internal::band_solve_triangular_selector<int,Upper|UnitDiag,Scalar,false,Scalar,ColMajor>::run);
352     func[TR    | (UP << 2) | (UNIT  << 3)] = (internal::band_solve_triangular_selector<int,Lower|UnitDiag,Scalar,false,Scalar,RowMajor>::run);
353     func[ADJ   | (UP << 2) | (UNIT  << 3)] = (internal::band_solve_triangular_selector<int,Lower|UnitDiag,Scalar,Conj, Scalar,RowMajor>::run);
354 
355     func[NOTR  | (LO << 2) | (UNIT  << 3)] = (internal::band_solve_triangular_selector<int,Lower|UnitDiag,Scalar,false,Scalar,ColMajor>::run);
356     func[TR    | (LO << 2) | (UNIT  << 3)] = (internal::band_solve_triangular_selector<int,Upper|UnitDiag,Scalar,false,Scalar,RowMajor>::run);
357     func[ADJ   | (LO << 2) | (UNIT  << 3)] = (internal::band_solve_triangular_selector<int,Upper|UnitDiag,Scalar,Conj, Scalar,RowMajor>::run);
358 
359     init = true;
360   }
361 
362   Scalar* a = reinterpret_cast<Scalar*>(pa);
363   Scalar* x = reinterpret_cast<Scalar*>(px);
364   int coeff_rows = *k+1;
365 
366   int info = 0;
367        if(UPLO(*uplo)==INVALID)                                       info = 1;
368   else if(OP(*op)==INVALID)                                           info = 2;
369   else if(DIAG(*diag)==INVALID)                                       info = 3;
370   else if(*n<0)                                                       info = 4;
371   else if(*k<0)                                                       info = 5;
372   else if(*lda<coeff_rows)                                            info = 7;
373   else if(*incx==0)                                                   info = 9;
374   if(info)
375     return xerbla_(SCALAR_SUFFIX_UP"TBSV ",&info,6);
376 
377   if(*n==0 || (*k==0 && DIAG(*diag)==UNIT))
378     return 0;
379 
380   int actual_n = *n;
381 
382   Scalar* actual_x = get_compact_vector(x,actual_n,*incx);
383 
384   int code = OP(*op) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
385   if(code>=16 || func[code]==0)
386     return 0;
387 
388   func[code](*n, *k, a, *lda, actual_x);
389 
390   if(actual_x!=x) delete[] copy_back(actual_x,x,actual_n,*incx);
391 
392   return 0;
393 }
394 
395 /**  DTPMV  performs one of the matrix-vector operations
396   *
397   *     x := A*x,   or   x := A'*x,
398   *
399   *  where x is an n element vector and  A is an n by n unit, or non-unit,
400   *  upper or lower triangular matrix, supplied in packed form.
401   */
EIGEN_BLAS_FUNC(tpmv)402 int EIGEN_BLAS_FUNC(tpmv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pap, RealScalar *px, int *incx)
403 {
404   typedef void (*functype)(int, const Scalar*, const Scalar*, Scalar*, Scalar);
405   static functype func[16];
406 
407   static bool init = false;
408   if(!init)
409   {
410     for(int k=0; k<16; ++k)
411       func[k] = 0;
412 
413     func[NOTR  | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|0,       Scalar,false,Scalar,false,ColMajor>::run);
414     func[TR    | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|0,       Scalar,false,Scalar,false,RowMajor>::run);
415     func[ADJ   | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|0,       Scalar,Conj, Scalar,false,RowMajor>::run);
416 
417     func[NOTR  | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|0,       Scalar,false,Scalar,false,ColMajor>::run);
418     func[TR    | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|0,       Scalar,false,Scalar,false,RowMajor>::run);
419     func[ADJ   | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|0,       Scalar,Conj, Scalar,false,RowMajor>::run);
420 
421     func[NOTR  | (UP << 2) | (UNIT  << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run);
422     func[TR    | (UP << 2) | (UNIT  << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run);
423     func[ADJ   | (UP << 2) | (UNIT  << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run);
424 
425     func[NOTR  | (LO << 2) | (UNIT  << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run);
426     func[TR    | (LO << 2) | (UNIT  << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run);
427     func[ADJ   | (LO << 2) | (UNIT  << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run);
428 
429     init = true;
430   }
431 
432   Scalar* ap = reinterpret_cast<Scalar*>(pap);
433   Scalar* x = reinterpret_cast<Scalar*>(px);
434 
435   int info = 0;
436   if(UPLO(*uplo)==INVALID)                                            info = 1;
437   else if(OP(*opa)==INVALID)                                          info = 2;
438   else if(DIAG(*diag)==INVALID)                                       info = 3;
439   else if(*n<0)                                                       info = 4;
440   else if(*incx==0)                                                   info = 7;
441   if(info)
442     return xerbla_(SCALAR_SUFFIX_UP"TPMV ",&info,6);
443 
444   if(*n==0)
445     return 1;
446 
447   Scalar* actual_x = get_compact_vector(x,*n,*incx);
448   Matrix<Scalar,Dynamic,1> res(*n);
449   res.setZero();
450 
451   int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
452   if(code>=16 || func[code]==0)
453     return 0;
454 
455   func[code](*n, ap, actual_x, res.data(), Scalar(1));
456 
457   copy_back(res.data(),x,*n,*incx);
458   if(actual_x!=x) delete[] actual_x;
459 
460   return 1;
461 }
462 
463 /**  DTPSV  solves one of the systems of equations
464   *
465   *     A*x = b,   or   A'*x = b,
466   *
467   *  where b and x are n element vectors and A is an n by n unit, or
468   *  non-unit, upper or lower triangular matrix, supplied in packed form.
469   *
470   *  No test for singularity or near-singularity is included in this
471   *  routine. Such tests must be performed before calling this routine.
472   */
EIGEN_BLAS_FUNC(tpsv)473 int EIGEN_BLAS_FUNC(tpsv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pap, RealScalar *px, int *incx)
474 {
475   typedef void (*functype)(int, const Scalar*, Scalar*);
476   static functype func[16];
477 
478   static bool init = false;
479   if(!init)
480   {
481     for(int k=0; k<16; ++k)
482       func[k] = 0;
483 
484     func[NOTR  | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0,       false,ColMajor>::run);
485     func[TR    | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0,       false,RowMajor>::run);
486     func[ADJ   | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0,       Conj, RowMajor>::run);
487 
488     func[NOTR  | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0,       false,ColMajor>::run);
489     func[TR    | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0,       false,RowMajor>::run);
490     func[ADJ   | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0,       Conj, RowMajor>::run);
491 
492     func[NOTR  | (UP << 2) | (UNIT  << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,ColMajor>::run);
493     func[TR    | (UP << 2) | (UNIT  << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,RowMajor>::run);
494     func[ADJ   | (UP << 2) | (UNIT  << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,Conj, RowMajor>::run);
495 
496     func[NOTR  | (LO << 2) | (UNIT  << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,ColMajor>::run);
497     func[TR    | (LO << 2) | (UNIT  << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,RowMajor>::run);
498     func[ADJ   | (LO << 2) | (UNIT  << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,Conj, RowMajor>::run);
499 
500     init = true;
501   }
502 
503   Scalar* ap = reinterpret_cast<Scalar*>(pap);
504   Scalar* x = reinterpret_cast<Scalar*>(px);
505 
506   int info = 0;
507   if(UPLO(*uplo)==INVALID)                                            info = 1;
508   else if(OP(*opa)==INVALID)                                          info = 2;
509   else if(DIAG(*diag)==INVALID)                                       info = 3;
510   else if(*n<0)                                                       info = 4;
511   else if(*incx==0)                                                   info = 7;
512   if(info)
513     return xerbla_(SCALAR_SUFFIX_UP"TPSV ",&info,6);
514 
515   Scalar* actual_x = get_compact_vector(x,*n,*incx);
516 
517   int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
518   func[code](*n, ap, actual_x);
519 
520   if(actual_x!=x) delete[] copy_back(actual_x,x,*n,*incx);
521 
522   return 1;
523 }
524 
525