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 }
SCALAR_SUFFIX_UP
#define SCALAR_SUFFIX_UP
Definition: blas/complex_double.cpp:12
m
Matrix3f m
Definition: AngleAxis_mimic_euler.cpp:1
general_matrix_vector_product_wrapper
Definition: level2_impl.h:13
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46
max
#define max(a, b)
Definition: datatypes.h:20
Eigen::PlainObjectBase::setZero
EIGEN_DEVICE_FUNC Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:562
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
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
y
Scalar * y
Definition: level1_cplx_impl.h:124
b
Scalar * b
Definition: benchVecAdd.cpp:17
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
gtsam::diag
Matrix diag(const std::vector< Matrix > &Hs)
Definition: Matrix.cpp:206
Eigen::PlainObjectBase::data
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar * data() const
Definition: PlainObjectBase.h:247
min
#define min(a, b)
Definition: datatypes.h:19
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:119
incx
RealScalar RealScalar int * incx
Definition: level1_cplx_impl.h:29
n
int n
Definition: BiCGSTAB_simple.cpp:1
c
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
copy_back
T * copy_back(T *x_cpy, T *x, int n, int incx)
Definition: gtsam/3rdparty/Eigen/blas/common.h:159
offset
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy y set format x g set format y g set format x2 g set format y2 g set format z g set angles radians set nogrid set key title set key left top Right noreverse box linetype linewidth samplen spacing width set nolabel set noarrow set nologscale set logscale x set set pointsize set encoding default set nopolar set noparametric set set set set surface set nocontour set clabel set mapping cartesian set nohidden3d set cntrparam order set cntrparam linear set cntrparam levels auto set cntrparam points set size set set xzeroaxis lt lw set x2zeroaxis lt lw set yzeroaxis lt lw set y2zeroaxis lt lw set tics in set ticslevel set tics set mxtics default set mytics default set mx2tics default set my2tics default set xtics border mirror norotate autofreq set ytics border mirror norotate autofreq set ztics border nomirror norotate autofreq set nox2tics set noy2tics set timestamp bottom norotate offset
Definition: gnuplot_common_settings.hh:12
MatrixType
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
trans
static char trans
Definition: blas_interface.hh:58
EIGEN_BLAS_FUNC
#define EIGEN_BLAS_FUNC(X)
Definition: gtsam/3rdparty/Eigen/blas/common.h:173
info
else if n * info
Definition: 3rdparty/Eigen/lapack/cholesky.cpp:18
res
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
Definition: PartialRedux_count.cpp:3
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
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
test_callbacks.value
value
Definition: test_callbacks.py:148
get_compact_vector
T * get_compact_vector(T *x, int n, int incx)
Definition: gtsam/3rdparty/Eigen/blas/common.h:147
py
int RealScalar int RealScalar * py
Definition: level1_cplx_impl.h:119
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
xerbla_
EIGEN_WEAK_LINKING int xerbla_(const char *msg, int *info, int)
Definition: xerbla.cpp:15
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
palpha
int RealScalar * palpha
Definition: level1_cplx_impl.h:142
internal::packed_triangular_solve_vector
Definition: PackedTriangularSolverVector.h:16
px
RealScalar RealScalar * px
Definition: level1_cplx_impl.h:28
alpha
RealScalar alpha
Definition: level1_cplx_impl.h:147
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
lda
* lda
Definition: eigenvalues.cpp:59
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
internal::band_solve_triangular_selector
Definition: BandTriangularSolver.h:19
std::swap
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)
Definition: NearestNeighbor.hpp:827
Eigen::placeholders::end
static EIGEN_DEPRECATED const end_t end
Definition: IndexedViewHelper.h:181
make_vector
Map< Matrix< T, Dynamic, 1 >, 0, InnerStride< Dynamic > > make_vector(T *data, int size, int incr)
Definition: gtsam/3rdparty/Eigen/blas/common.h:123
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
a
ArrayXXi a
Definition: Array_initializer_list_23_cxx11.cpp:1
len
size_t len(handle h)
Get the length of a Python object.
Definition: pytypes.h:2244
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2
incy
int RealScalar int RealScalar int * incy
Definition: level1_cplx_impl.h:119
internal::packed_triangular_matrix_vector_product
Definition: PackedTriangularMatrixVector.h:16
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
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

gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:32