level2_impl.h
Go to the documentation of this file.
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 
12 template<typename Index, typename Scalar, int StorageOrder, bool ConjugateLhs, bool ConjugateRhs>
13 struct general_matrix_vector_product_wrapper
14 {
15  static void run(Index rows, Index cols,const Scalar *lhs, Index lhsStride, const Scalar *rhs, Index rhsIncr, Scalar* res, Index resIncr, Scalar alpha)
16  {
17  typedef internal::const_blas_data_mapper<Scalar,Index,StorageOrder> LhsMapper;
18  typedef internal::const_blas_data_mapper<Scalar,Index,RowMajor> RhsMapper;
19 
20  internal::general_matrix_vector_product
21  <Index,Scalar,LhsMapper,StorageOrder,ConjugateLhs,Scalar,RhsMapper,ConjugateRhs>::run(
22  rows, cols, LhsMapper(lhs, lhsStride), RhsMapper(rhs, rhsIncr), res, resIncr, alpha);
23  }
24 };
25 
26 int EIGEN_BLAS_FUNC(gemv)(const char *opa, const int *m, const int *n, const RealScalar *palpha,
27  const RealScalar *pa, const int *lda, const RealScalar *pb, const int *incb, const RealScalar *pbeta, RealScalar *pc, const int *incc)
28 {
29  typedef void (*functype)(int, int, const Scalar *, int, const Scalar *, int , Scalar *, int, Scalar);
30  static const functype func[4] = {
31  // array index: NOTR
32  (general_matrix_vector_product_wrapper<int,Scalar,ColMajor,false,false>::run),
33  // array index: TR
34  (general_matrix_vector_product_wrapper<int,Scalar,RowMajor,false,false>::run),
35  // array index: ADJ
36  (general_matrix_vector_product_wrapper<int,Scalar,RowMajor,Conj ,false>::run),
37  0
38  };
39 
40  const Scalar* a = reinterpret_cast<const Scalar*>(pa);
41  const Scalar* b = reinterpret_cast<const Scalar*>(pb);
42  Scalar* c = reinterpret_cast<Scalar*>(pc);
43  Scalar alpha = *reinterpret_cast<const Scalar*>(palpha);
44  Scalar beta = *reinterpret_cast<const Scalar*>(pbeta);
45 
46  // check arguments
47  int info = 0;
48  if(OP(*opa)==INVALID) info = 1;
49  else if(*m<0) info = 2;
50  else if(*n<0) info = 3;
51  else if(*lda<std::max(1,*m)) info = 6;
52  else if(*incb==0) info = 8;
53  else if(*incc==0) info = 11;
54  if(info)
55  return xerbla_(SCALAR_SUFFIX_UP"GEMV ",&info,6);
56 
57  if(*m==0 || *n==0 || (alpha==Scalar(0) && beta==Scalar(1)))
58  return 0;
59 
60  int actual_m = *m;
61  int actual_n = *n;
62  int code = OP(*opa);
63  if(code!=NOTR)
64  std::swap(actual_m,actual_n);
65 
66  const Scalar* actual_b = get_compact_vector(b,actual_n,*incb);
67  Scalar* actual_c = get_compact_vector(c,actual_m,*incc);
68 
69  if(beta!=Scalar(1))
70  {
71  if(beta==Scalar(0)) make_vector(actual_c, actual_m).setZero();
72  else make_vector(actual_c, actual_m) *= beta;
73  }
74 
75  if(code>=4 || func[code]==0)
76  return 0;
77 
78  func[code](actual_m, actual_n, a, *lda, actual_b, 1, actual_c, 1, alpha);
79 
80  if(actual_b!=b) delete[] actual_b;
81  if(actual_c!=c) delete[] copy_back(actual_c,c,actual_m,*incc);
82 
83  return 1;
84 }
85 
86 int EIGEN_BLAS_FUNC(trsv)(const char *uplo, const char *opa, const char *diag, const int *n, const RealScalar *pa, const int *lda, RealScalar *pb, const int *incb)
87 {
88  typedef void (*functype)(int, const Scalar *, int, Scalar *);
89  static const functype func[16] = {
90  // array index: NOTR | (UP << 2) | (NUNIT << 3)
91  (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0, false,ColMajor>::run),
92  // array index: TR | (UP << 2) | (NUNIT << 3)
93  (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0, false,RowMajor>::run),
94  // array index: ADJ | (UP << 2) | (NUNIT << 3)
95  (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0, Conj, RowMajor>::run),
96  0,
97  // array index: NOTR | (LO << 2) | (NUNIT << 3)
98  (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0, false,ColMajor>::run),
99  // array index: TR | (LO << 2) | (NUNIT << 3)
100  (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0, false,RowMajor>::run),
101  // array index: ADJ | (LO << 2) | (NUNIT << 3)
102  (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0, Conj, RowMajor>::run),
103  0,
104  // array index: NOTR | (UP << 2) | (UNIT << 3)
105  (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,ColMajor>::run),
106  // array index: TR | (UP << 2) | (UNIT << 3)
107  (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,RowMajor>::run),
108  // array index: ADJ | (UP << 2) | (UNIT << 3)
109  (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,Conj, RowMajor>::run),
110  0,
111  // array index: NOTR | (LO << 2) | (UNIT << 3)
112  (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,ColMajor>::run),
113  // array index: TR | (LO << 2) | (UNIT << 3)
114  (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,RowMajor>::run),
115  // array index: ADJ | (LO << 2) | (UNIT << 3)
116  (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,Conj, RowMajor>::run),
117  0
118  };
119 
120  const Scalar* a = reinterpret_cast<const Scalar*>(pa);
121  Scalar* b = reinterpret_cast<Scalar*>(pb);
122 
123  int info = 0;
124  if(UPLO(*uplo)==INVALID) info = 1;
125  else if(OP(*opa)==INVALID) info = 2;
126  else if(DIAG(*diag)==INVALID) info = 3;
127  else if(*n<0) info = 4;
128  else if(*lda<std::max(1,*n)) info = 6;
129  else if(*incb==0) info = 8;
130  if(info)
131  return xerbla_(SCALAR_SUFFIX_UP"TRSV ",&info,6);
132 
133  Scalar* actual_b = get_compact_vector(b,*n,*incb);
134 
135  int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
136  func[code](*n, a, *lda, actual_b);
137 
138  if(actual_b!=b) delete[] copy_back(actual_b,b,*n,*incb);
139 
140  return 0;
141 }
142 
143 
144 
145 int EIGEN_BLAS_FUNC(trmv)(const char *uplo, const char *opa, const char *diag, const int *n, const RealScalar *pa, const int *lda, RealScalar *pb, const int *incb)
146 {
147  typedef void (*functype)(int, int, const Scalar *, int, const Scalar *, int, Scalar *, int, const Scalar&);
148  static const functype func[16] = {
149  // array index: NOTR | (UP << 2) | (NUNIT << 3)
150  (internal::triangular_matrix_vector_product<int,Upper|0, Scalar,false,Scalar,false,ColMajor>::run),
151  // array index: TR | (UP << 2) | (NUNIT << 3)
152  (internal::triangular_matrix_vector_product<int,Lower|0, Scalar,false,Scalar,false,RowMajor>::run),
153  // array index: ADJ | (UP << 2) | (NUNIT << 3)
154  (internal::triangular_matrix_vector_product<int,Lower|0, Scalar,Conj, Scalar,false,RowMajor>::run),
155  0,
156  // array index: NOTR | (LO << 2) | (NUNIT << 3)
157  (internal::triangular_matrix_vector_product<int,Lower|0, Scalar,false,Scalar,false,ColMajor>::run),
158  // array index: TR | (LO << 2) | (NUNIT << 3)
159  (internal::triangular_matrix_vector_product<int,Upper|0, Scalar,false,Scalar,false,RowMajor>::run),
160  // array index: ADJ | (LO << 2) | (NUNIT << 3)
161  (internal::triangular_matrix_vector_product<int,Upper|0, Scalar,Conj, Scalar,false,RowMajor>::run),
162  0,
163  // array index: NOTR | (UP << 2) | (UNIT << 3)
164  (internal::triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run),
165  // array index: TR | (UP << 2) | (UNIT << 3)
166  (internal::triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run),
167  // array index: ADJ | (UP << 2) | (UNIT << 3)
168  (internal::triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run),
169  0,
170  // array index: NOTR | (LO << 2) | (UNIT << 3)
171  (internal::triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run),
172  // array index: TR | (LO << 2) | (UNIT << 3)
173  (internal::triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run),
174  // array index: ADJ | (LO << 2) | (UNIT << 3)
175  (internal::triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run),
176  0
177  };
178 
179  const Scalar* a = reinterpret_cast<const Scalar*>(pa);
180  Scalar* b = reinterpret_cast<Scalar*>(pb);
181 
182  int info = 0;
183  if(UPLO(*uplo)==INVALID) info = 1;
184  else if(OP(*opa)==INVALID) info = 2;
185  else if(DIAG(*diag)==INVALID) info = 3;
186  else if(*n<0) info = 4;
187  else if(*lda<std::max(1,*n)) info = 6;
188  else if(*incb==0) info = 8;
189  if(info)
190  return xerbla_(SCALAR_SUFFIX_UP"TRMV ",&info,6);
191 
192  if(*n==0)
193  return 1;
194 
195  Scalar* actual_b = get_compact_vector(b,*n,*incb);
196  Matrix<Scalar,Dynamic,1> res(*n);
197  res.setZero();
198 
199  int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
200  if(code>=16 || func[code]==0)
201  return 0;
202 
203  func[code](*n, *n, a, *lda, actual_b, 1, res.data(), 1, Scalar(1));
204 
205  copy_back(res.data(),b,*n,*incb);
206  if(actual_b!=b) delete[] actual_b;
207 
208  return 1;
209 }
210 
218 int EIGEN_BLAS_FUNC(gbmv)(char *trans, int *m, int *n, int *kl, int *ku, RealScalar *palpha, RealScalar *pa, int *lda,
219  RealScalar *px, int *incx, RealScalar *pbeta, RealScalar *py, int *incy)
220 {
221  const Scalar* a = reinterpret_cast<const Scalar*>(pa);
222  const Scalar* x = reinterpret_cast<const Scalar*>(px);
223  Scalar* y = reinterpret_cast<Scalar*>(py);
224  Scalar alpha = *reinterpret_cast<const Scalar*>(palpha);
225  Scalar beta = *reinterpret_cast<const Scalar*>(pbeta);
226  int coeff_rows = *kl+*ku+1;
227 
228  int info = 0;
229  if(OP(*trans)==INVALID) info = 1;
230  else if(*m<0) info = 2;
231  else if(*n<0) info = 3;
232  else if(*kl<0) info = 4;
233  else if(*ku<0) info = 5;
234  else if(*lda<coeff_rows) info = 8;
235  else if(*incx==0) info = 10;
236  else if(*incy==0) info = 13;
237  if(info)
238  return xerbla_(SCALAR_SUFFIX_UP"GBMV ",&info,6);
239 
240  if(*m==0 || *n==0 || (alpha==Scalar(0) && beta==Scalar(1)))
241  return 0;
242 
243  int actual_m = *m;
244  int actual_n = *n;
245  if(OP(*trans)!=NOTR)
246  std::swap(actual_m,actual_n);
247 
248  const Scalar* actual_x = get_compact_vector(x,actual_n,*incx);
249  Scalar* actual_y = get_compact_vector(y,actual_m,*incy);
250 
251  if(beta!=Scalar(1))
252  {
253  if(beta==Scalar(0)) make_vector(actual_y, actual_m).setZero();
254  else make_vector(actual_y, actual_m) *= beta;
255  }
256 
257  ConstMatrixType mat_coeffs(a,coeff_rows,*n,*lda);
258 
259  int nb = std::min(*n,(*m)+(*ku));
260  for(int j=0; j<nb; ++j)
261  {
262  int start = std::max(0,j - *ku);
263  int end = std::min((*m)-1,j + *kl);
264  int len = end - start + 1;
265  int offset = (*ku) - j + start;
266  if(OP(*trans)==NOTR)
267  make_vector(actual_y+start,len) += (alpha*actual_x[j]) * mat_coeffs.col(j).segment(offset,len);
268  else if(OP(*trans)==TR)
269  actual_y[j] += alpha * ( mat_coeffs.col(j).segment(offset,len).transpose() * make_vector(actual_x+start,len) ).value();
270  else
271  actual_y[j] += alpha * ( mat_coeffs.col(j).segment(offset,len).adjoint() * make_vector(actual_x+start,len) ).value();
272  }
273 
274  if(actual_x!=x) delete[] actual_x;
275  if(actual_y!=y) delete[] copy_back(actual_y,y,actual_m,*incy);
276 
277  return 0;
278 }
279 
280 #if 0
281 
288 int EIGEN_BLAS_FUNC(tbmv)(char *uplo, char *opa, char *diag, int *n, int *k, RealScalar *pa, int *lda, RealScalar *px, int *incx)
289 {
290  Scalar* a = reinterpret_cast<Scalar*>(pa);
291  Scalar* x = reinterpret_cast<Scalar*>(px);
292  int coeff_rows = *k + 1;
293 
294  int info = 0;
295  if(UPLO(*uplo)==INVALID) info = 1;
296  else if(OP(*opa)==INVALID) info = 2;
297  else if(DIAG(*diag)==INVALID) info = 3;
298  else if(*n<0) info = 4;
299  else if(*k<0) info = 5;
300  else if(*lda<coeff_rows) info = 7;
301  else if(*incx==0) info = 9;
302  if(info)
303  return xerbla_(SCALAR_SUFFIX_UP"TBMV ",&info,6);
304 
305  if(*n==0)
306  return 0;
307 
308  int actual_n = *n;
309 
310  Scalar* actual_x = get_compact_vector(x,actual_n,*incx);
311 
312  MatrixType mat_coeffs(a,coeff_rows,*n,*lda);
313 
314  int ku = UPLO(*uplo)==UPPER ? *k : 0;
315  int kl = UPLO(*uplo)==LOWER ? *k : 0;
316 
317  for(int j=0; j<*n; ++j)
318  {
319  int start = std::max(0,j - ku);
320  int end = std::min((*m)-1,j + kl);
321  int len = end - start + 1;
322  int offset = (ku) - j + start;
323 
324  if(OP(*trans)==NOTR)
325  make_vector(actual_y+start,len) += (alpha*actual_x[j]) * mat_coeffs.col(j).segment(offset,len);
326  else if(OP(*trans)==TR)
327  actual_y[j] += alpha * ( mat_coeffs.col(j).segment(offset,len).transpose() * make_vector(actual_x+start,len) ).value();
328  else
329  actual_y[j] += alpha * ( mat_coeffs.col(j).segment(offset,len).adjoint() * make_vector(actual_x+start,len) ).value();
330  }
331 
332  if(actual_x!=x) delete[] actual_x;
333  if(actual_y!=y) delete[] copy_back(actual_y,y,actual_m,*incy);
334 
335  return 0;
336 }
337 #endif
338 
350 int EIGEN_BLAS_FUNC(tbsv)(char *uplo, char *op, char *diag, int *n, int *k, RealScalar *pa, int *lda, RealScalar *px, int *incx)
351 {
352  typedef void (*functype)(int, int, const Scalar *, int, Scalar *);
353  static const functype func[16] = {
354  // array index: NOTR | (UP << 2) | (NUNIT << 3)
355  (internal::band_solve_triangular_selector<int,Upper|0, Scalar,false,Scalar,ColMajor>::run),
356  // array index: TR | (UP << 2) | (NUNIT << 3)
357  (internal::band_solve_triangular_selector<int,Lower|0, Scalar,false,Scalar,RowMajor>::run),
358  // array index: ADJ | (UP << 2) | (NUNIT << 3)
359  (internal::band_solve_triangular_selector<int,Lower|0, Scalar,Conj, Scalar,RowMajor>::run),
360  0,
361  // array index: NOTR | (LO << 2) | (NUNIT << 3)
362  (internal::band_solve_triangular_selector<int,Lower|0, Scalar,false,Scalar,ColMajor>::run),
363  // array index: TR | (LO << 2) | (NUNIT << 3)
364  (internal::band_solve_triangular_selector<int,Upper|0, Scalar,false,Scalar,RowMajor>::run),
365  // array index: ADJ | (LO << 2) | (NUNIT << 3)
366  (internal::band_solve_triangular_selector<int,Upper|0, Scalar,Conj, Scalar,RowMajor>::run),
367  0,
368  // array index: NOTR | (UP << 2) | (UNIT << 3)
369  (internal::band_solve_triangular_selector<int,Upper|UnitDiag,Scalar,false,Scalar,ColMajor>::run),
370  // array index: TR | (UP << 2) | (UNIT << 3)
371  (internal::band_solve_triangular_selector<int,Lower|UnitDiag,Scalar,false,Scalar,RowMajor>::run),
372  // array index: ADJ | (UP << 2) | (UNIT << 3)
373  (internal::band_solve_triangular_selector<int,Lower|UnitDiag,Scalar,Conj, Scalar,RowMajor>::run),
374  0,
375  // array index: NOTR | (LO << 2) | (UNIT << 3)
376  (internal::band_solve_triangular_selector<int,Lower|UnitDiag,Scalar,false,Scalar,ColMajor>::run),
377  // array index: TR | (LO << 2) | (UNIT << 3)
378  (internal::band_solve_triangular_selector<int,Upper|UnitDiag,Scalar,false,Scalar,RowMajor>::run),
379  // array index: ADJ | (LO << 2) | (UNIT << 3)
380  (internal::band_solve_triangular_selector<int,Upper|UnitDiag,Scalar,Conj, Scalar,RowMajor>::run),
381  0,
382  };
383 
384  Scalar* a = reinterpret_cast<Scalar*>(pa);
385  Scalar* x = reinterpret_cast<Scalar*>(px);
386  int coeff_rows = *k+1;
387 
388  int info = 0;
389  if(UPLO(*uplo)==INVALID) info = 1;
390  else if(OP(*op)==INVALID) info = 2;
391  else if(DIAG(*diag)==INVALID) info = 3;
392  else if(*n<0) info = 4;
393  else if(*k<0) info = 5;
394  else if(*lda<coeff_rows) info = 7;
395  else if(*incx==0) info = 9;
396  if(info)
397  return xerbla_(SCALAR_SUFFIX_UP"TBSV ",&info,6);
398 
399  if(*n==0 || (*k==0 && DIAG(*diag)==UNIT))
400  return 0;
401 
402  int actual_n = *n;
403 
404  Scalar* actual_x = get_compact_vector(x,actual_n,*incx);
405 
406  int code = OP(*op) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
407  if(code>=16 || func[code]==0)
408  return 0;
409 
410  func[code](*n, *k, a, *lda, actual_x);
411 
412  if(actual_x!=x) delete[] copy_back(actual_x,x,actual_n,*incx);
413 
414  return 0;
415 }
416 
424 int EIGEN_BLAS_FUNC(tpmv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pap, RealScalar *px, int *incx)
425 {
426  typedef void (*functype)(int, const Scalar*, const Scalar*, Scalar*, Scalar);
427  static const functype func[16] = {
428  // array index: NOTR | (UP << 2) | (NUNIT << 3)
429  (internal::packed_triangular_matrix_vector_product<int,Upper|0, Scalar,false,Scalar,false,ColMajor>::run),
430  // array index: TR | (UP << 2) | (NUNIT << 3)
431  (internal::packed_triangular_matrix_vector_product<int,Lower|0, Scalar,false,Scalar,false,RowMajor>::run),
432  // array index: ADJ | (UP << 2) | (NUNIT << 3)
433  (internal::packed_triangular_matrix_vector_product<int,Lower|0, Scalar,Conj, Scalar,false,RowMajor>::run),
434  0,
435  // array index: NOTR | (LO << 2) | (NUNIT << 3)
436  (internal::packed_triangular_matrix_vector_product<int,Lower|0, Scalar,false,Scalar,false,ColMajor>::run),
437  // array index: TR | (LO << 2) | (NUNIT << 3)
438  (internal::packed_triangular_matrix_vector_product<int,Upper|0, Scalar,false,Scalar,false,RowMajor>::run),
439  // array index: ADJ | (LO << 2) | (NUNIT << 3)
440  (internal::packed_triangular_matrix_vector_product<int,Upper|0, Scalar,Conj, Scalar,false,RowMajor>::run),
441  0,
442  // array index: NOTR | (UP << 2) | (UNIT << 3)
443  (internal::packed_triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run),
444  // array index: TR | (UP << 2) | (UNIT << 3)
445  (internal::packed_triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run),
446  // array index: ADJ | (UP << 2) | (UNIT << 3)
447  (internal::packed_triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run),
448  0,
449  // array index: NOTR | (LO << 2) | (UNIT << 3)
450  (internal::packed_triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run),
451  // array index: TR | (LO << 2) | (UNIT << 3)
452  (internal::packed_triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run),
453  // array index: ADJ | (LO << 2) | (UNIT << 3)
454  (internal::packed_triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run),
455  0
456  };
457 
458  Scalar* ap = reinterpret_cast<Scalar*>(pap);
459  Scalar* x = reinterpret_cast<Scalar*>(px);
460 
461  int info = 0;
462  if(UPLO(*uplo)==INVALID) info = 1;
463  else if(OP(*opa)==INVALID) info = 2;
464  else if(DIAG(*diag)==INVALID) info = 3;
465  else if(*n<0) info = 4;
466  else if(*incx==0) info = 7;
467  if(info)
468  return xerbla_(SCALAR_SUFFIX_UP"TPMV ",&info,6);
469 
470  if(*n==0)
471  return 1;
472 
473  Scalar* actual_x = get_compact_vector(x,*n,*incx);
474  Matrix<Scalar,Dynamic,1> res(*n);
475  res.setZero();
476 
477  int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
478  if(code>=16 || func[code]==0)
479  return 0;
480 
481  func[code](*n, ap, actual_x, res.data(), Scalar(1));
482 
483  copy_back(res.data(),x,*n,*incx);
484  if(actual_x!=x) delete[] actual_x;
485 
486  return 1;
487 }
488 
499 int EIGEN_BLAS_FUNC(tpsv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pap, RealScalar *px, int *incx)
500 {
501  typedef void (*functype)(int, const Scalar*, Scalar*);
502  static const functype func[16] = {
503  // array index: NOTR | (UP << 2) | (NUNIT << 3)
504  (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0, false,ColMajor>::run),
505  // array index: TR | (UP << 2) | (NUNIT << 3)
506  (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0, false,RowMajor>::run),
507  // array index: ADJ | (UP << 2) | (NUNIT << 3)
508  (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0, Conj, RowMajor>::run),
509  0,
510  // array index: NOTR | (LO << 2) | (NUNIT << 3)
511  (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0, false,ColMajor>::run),
512  // array index: TR | (LO << 2) | (NUNIT << 3)
513  (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0, false,RowMajor>::run),
514  // array index: ADJ | (LO << 2) | (NUNIT << 3)
515  (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0, Conj, RowMajor>::run),
516  0,
517  // array index: NOTR | (UP << 2) | (UNIT << 3)
518  (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,ColMajor>::run),
519  // array index: TR | (UP << 2) | (UNIT << 3)
520  (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,RowMajor>::run),
521  // array index: ADJ | (UP << 2) | (UNIT << 3)
522  (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,Conj, RowMajor>::run),
523  0,
524  // array index: NOTR | (LO << 2) | (UNIT << 3)
525  (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,ColMajor>::run),
526  // array index: TR | (LO << 2) | (UNIT << 3)
527  (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,RowMajor>::run),
528  // array index: ADJ | (LO << 2) | (UNIT << 3)
529  (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,Conj, RowMajor>::run),
530  0
531  };
532 
533  Scalar* ap = reinterpret_cast<Scalar*>(pap);
534  Scalar* x = reinterpret_cast<Scalar*>(px);
535 
536  int info = 0;
537  if(UPLO(*uplo)==INVALID) info = 1;
538  else if(OP(*opa)==INVALID) info = 2;
539  else if(DIAG(*diag)==INVALID) info = 3;
540  else if(*n<0) info = 4;
541  else if(*incx==0) info = 7;
542  if(info)
543  return xerbla_(SCALAR_SUFFIX_UP"TPSV ",&info,6);
544 
545  Scalar* actual_x = get_compact_vector(x,*n,*incx);
546 
547  int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
548  func[code](*n, ap, actual_x);
549 
550  if(actual_x!=x) delete[] copy_back(actual_x,x,*n,*incx);
551 
552  return 1;
553 }
trmv
int EIGEN_BLAS_FUNC() trmv(const char *uplo, const char *opa, const char *diag, const int *n, const RealScalar *pa, const int *lda, RealScalar *pb, const int *incb)
Definition: level2_impl.h:145
DIAG
#define DIAG(X)
Definition: common.h:52
UNIT
#define UNIT
Definition: common.h:35
b
Scalar * b
Definition: cholesky.cpp:56
alpha
RealScalar alpha
Definition: level1_cplx_impl.h:125
gemv
int EIGEN_BLAS_FUNC() gemv(const char *opa, const int *m, const int *n, const RealScalar *palpha, const RealScalar *pa, const int *lda, const RealScalar *pb, const int *incb, const RealScalar *pbeta, RealScalar *pc, const int *incc)
Definition: level2_impl.h:26
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: common.h:85
TR
#define TR
Definition: common.h:25
internal::band_solve_triangular_selector
Definition: BandTriangularSolver.h:19
tpmv
int EIGEN_BLAS_FUNC() tpmv(char *uplo, char *opa, char *diag, int *n, RealScalar *pap, RealScalar *px, int *incx)
Definition: level2_impl.h:424
copy_back
T * copy_back(T *x_cpy, T *x, int n, int incx)
Definition: common.h:151
Scalar
SCALAR Scalar
Definition: common.h:84
get_compact_vector
T * get_compact_vector(T *x, int n, int incx)
Definition: common.h:139
INVALID
#define INVALID
Definition: common.h:37
EIGEN_BLAS_FUNC
#define EIGEN_BLAS_FUNC(X)
Definition: common.h:161
make_vector
Map< Matrix< T, Dynamic, 1 >, 0, InnerStride< Dynamic > > make_vector(T *data, int size, int incr)
Definition: common.h:115
info
else if n * info
Definition: cholesky.cpp:18
general_matrix_vector_product_wrapper
Definition: level2_impl.h:13
incy
int RealScalar int RealScalar int * incy
Definition: level1_cplx_impl.h:97
general_matrix_vector_product_wrapper::run
static void run(Index rows, Index cols, const Scalar *lhs, Index lhsStride, const Scalar *rhs, Index rhsIncr, Scalar *res, Index resIncr, Scalar alpha)
Definition: level2_impl.h:15
x
Scalar * x
Definition: level1_cplx_impl.h:89
incx
RealScalar RealScalar int * incx
Definition: level1_cplx_impl.h:29
SCALAR_SUFFIX_UP
#define SCALAR_SUFFIX_UP
Definition: blas/complex_double.cpp:12
internal::packed_triangular_matrix_vector_product
Definition: PackedTriangularMatrixVector.h:16
Eigen::Map
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
pc
int RealScalar int RealScalar int RealScalar * pc
Definition: level1_cplx_impl.h:97
lda
* lda
Definition: eigenvalues.cpp:59
int
return int(ret)+1
y
Scalar * y
Definition: level1_cplx_impl.h:102
OP
#define OP(X)
Definition: common.h:39
xerbla_
EIGEN_WEAK_LINKING int xerbla_(const char *msg, int *info, int)
Definition: xerbla.cpp:15
c
RealScalar c
Definition: level1_cplx_impl.h:103
common.h
a
Scalar * a
Definition: cholesky.cpp:26
Eigen::TensorSycl::run
void run(Expr &expr, Dev &dev)
Definition: TensorSyclRun.h:47
trsv
int EIGEN_BLAS_FUNC() trsv(const char *uplo, const char *opa, const char *diag, const int *n, const RealScalar *pa, const int *lda, RealScalar *pb, const int *incb)
Definition: level2_impl.h:86
Eigen::PlainObjectBase::data
EIGEN_DEVICE_FUNC const EIGEN_STRONG_INLINE Scalar * data() const
Definition: PlainObjectBase.h:255
gbmv
int EIGEN_BLAS_FUNC() gbmv(char *trans, int *m, int *n, int *kl, int *ku, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *px, int *incx, RealScalar *pbeta, RealScalar *py, int *incy)
Definition: level2_impl.h:218
min
#define min(a, b)
Definition: datatypes.h:19
tpsv
int EIGEN_BLAS_FUNC() tpsv(char *uplo, char *opa, char *diag, int *n, RealScalar *pap, RealScalar *px, int *incx)
Definition: level2_impl.h:499
Eigen::Matrix< Scalar, Dynamic, 1 >
py
int RealScalar int RealScalar * py
Definition: level1_cplx_impl.h:97
Eigen::PlainObjectBase::setZero
EIGEN_DEVICE_FUNC Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:515
n
PlainMatrixType mat * n
Definition: eigenvalues.cpp:41
palpha
int RealScalar * palpha
Definition: level1_cplx_impl.h:120
UPLO
#define UPLO(X)
Definition: common.h:48
max
#define max(a, b)
Definition: datatypes.h:20
swap
int EIGEN_BLAS_FUNC() swap(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
Definition: level1_impl.h:152
NOTR
#define NOTR
Definition: common.h:24
internal::packed_triangular_solve_vector
Definition: PackedTriangularSolverVector.h:16
tbsv
int EIGEN_BLAS_FUNC() tbsv(char *uplo, char *op, char *diag, int *n, int *k, RealScalar *pa, int *lda, RealScalar *px, int *incx)
Definition: level2_impl.h:350
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
px
RealScalar RealScalar * px
Definition: level1_cplx_impl.h:28

control_box_rst
Author(s): Christoph Rösmann
autogenerated on Wed Mar 2 2022 00:05:52