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  static void run(Index rows, Index cols, const Scalar *lhs, Index lhsStride, const Scalar *rhs, Index rhsIncr,
15  Scalar *res, Index resIncr, Scalar alpha) {
16  typedef Eigen::internal::const_blas_data_mapper<Scalar, Index, StorageOrder> LhsMapper;
17  typedef Eigen::internal::const_blas_data_mapper<Scalar, Index, Eigen::RowMajor> RhsMapper;
18 
19  Eigen::internal::general_matrix_vector_product<Index, Scalar, LhsMapper, StorageOrder, ConjugateLhs, Scalar,
20  RhsMapper, ConjugateRhs>::run(rows, cols, LhsMapper(lhs, lhsStride),
21  RhsMapper(rhs, rhsIncr), res, resIncr,
22  alpha);
23  }
24 };
25 
26 EIGEN_BLAS_FUNC(gemv)
27 (const char *opa, const int *m, const int *n, const RealScalar *palpha, const RealScalar *pa, const int *lda,
28  const RealScalar *pb, const int *incb, const RealScalar *pbeta, RealScalar *pc, const int *incc) {
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, Eigen::ColMajor, false, false>::run),
33  // array index: TR
34  (general_matrix_vector_product_wrapper<int, Scalar, Eigen::RowMajor, false, false>::run),
35  // array index: ADJ
36  (general_matrix_vector_product_wrapper<int, Scalar, Eigen::RowMajor, Conj, false>::run), 0};
37 
38  const Scalar *a = reinterpret_cast<const Scalar *>(pa);
39  const Scalar *b = reinterpret_cast<const Scalar *>(pb);
40  Scalar *c = reinterpret_cast<Scalar *>(pc);
41  Scalar alpha = *reinterpret_cast<const Scalar *>(palpha);
42  Scalar beta = *reinterpret_cast<const Scalar *>(pbeta);
43 
44  // check arguments
45  int info = 0;
46  if (OP(*opa) == INVALID)
47  info = 1;
48  else if (*m < 0)
49  info = 2;
50  else if (*n < 0)
51  info = 3;
52  else if (*lda < std::max(1, *m))
53  info = 6;
54  else if (*incb == 0)
55  info = 8;
56  else if (*incc == 0)
57  info = 11;
58  if (info) return xerbla_(SCALAR_SUFFIX_UP "GEMV ", &info);
59 
60  if (*m == 0 || *n == 0 || (alpha == Scalar(0) && beta == Scalar(1))) return;
61 
62  int actual_m = *m;
63  int actual_n = *n;
64  int code = OP(*opa);
65  if (code != NOTR) std::swap(actual_m, actual_n);
66 
67  const Scalar *actual_b = get_compact_vector(b, actual_n, *incb);
68  Scalar *actual_c = get_compact_vector(c, actual_m, *incc);
69 
70  if (beta != Scalar(1)) {
71  if (beta == Scalar(0))
72  make_vector(actual_c, actual_m).setZero();
73  else
74  make_vector(actual_c, actual_m) *= beta;
75  }
76 
77  if (code >= 4 || func[code] == 0) return;
78 
79  func[code](actual_m, actual_n, a, *lda, actual_b, 1, actual_c, 1, alpha);
80 
81  if (actual_b != b) delete[] actual_b;
82  if (actual_c != c) delete[] copy_back(actual_c, c, actual_m, *incc);
83 }
84 
85 EIGEN_BLAS_FUNC(trsv)
86 (const char *uplo, const char *opa, const char *diag, const int *n, const RealScalar *pa, const int *lda,
87  RealScalar *pb, const int *incb) {
88  typedef void (*functype)(int, const Scalar *, int, Scalar *);
89  using Eigen::ColMajor;
90  using Eigen::Lower;
91  using Eigen::OnTheLeft;
92  using Eigen::RowMajor;
93  using Eigen::UnitDiag;
94  using Eigen::Upper;
95  static const functype func[16] = {
96  // array index: NOTR | (UP << 2) | (NUNIT << 3)
97  (Eigen::internal::triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Upper | 0, false, ColMajor>::run),
98  // array index: TR | (UP << 2) | (NUNIT << 3)
99  (Eigen::internal::triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Lower | 0, false, RowMajor>::run),
100  // array index: ADJ | (UP << 2) | (NUNIT << 3)
101  (Eigen::internal::triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Lower | 0, Conj, RowMajor>::run), 0,
102  // array index: NOTR | (LO << 2) | (NUNIT << 3)
103  (Eigen::internal::triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Lower | 0, false, ColMajor>::run),
104  // array index: TR | (LO << 2) | (NUNIT << 3)
105  (Eigen::internal::triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Upper | 0, false, RowMajor>::run),
106  // array index: ADJ | (LO << 2) | (NUNIT << 3)
107  (Eigen::internal::triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Upper | 0, Conj, RowMajor>::run), 0,
108  // array index: NOTR | (UP << 2) | (UNIT << 3)
109  (Eigen::internal::triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Upper | UnitDiag, false,
110  ColMajor>::run),
111  // array index: TR | (UP << 2) | (UNIT << 3)
112  (Eigen::internal::triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Lower | UnitDiag, false,
113  RowMajor>::run),
114  // array index: ADJ | (UP << 2) | (UNIT << 3)
115  (Eigen::internal::triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Lower | UnitDiag, Conj, RowMajor>::run),
116  0,
117  // array index: NOTR | (LO << 2) | (UNIT << 3)
118  (Eigen::internal::triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Lower | UnitDiag, false,
119  ColMajor>::run),
120  // array index: TR | (LO << 2) | (UNIT << 3)
121  (Eigen::internal::triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Upper | UnitDiag, false,
122  RowMajor>::run),
123  // array index: ADJ | (LO << 2) | (UNIT << 3)
124  (Eigen::internal::triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Upper | UnitDiag, Conj, RowMajor>::run),
125  0};
126 
127  const Scalar *a = reinterpret_cast<const Scalar *>(pa);
128  Scalar *b = reinterpret_cast<Scalar *>(pb);
129 
130  int info = 0;
131  if (UPLO(*uplo) == INVALID)
132  info = 1;
133  else if (OP(*opa) == INVALID)
134  info = 2;
135  else if (DIAG(*diag) == INVALID)
136  info = 3;
137  else if (*n < 0)
138  info = 4;
139  else if (*lda < std::max(1, *n))
140  info = 6;
141  else if (*incb == 0)
142  info = 8;
143  if (info) return xerbla_(SCALAR_SUFFIX_UP "TRSV ", &info);
144 
145  Scalar *actual_b = get_compact_vector(b, *n, *incb);
146 
147  int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
148  func[code](*n, a, *lda, actual_b);
149 
150  if (actual_b != b) delete[] copy_back(actual_b, b, *n, *incb);
151 }
152 
153 EIGEN_BLAS_FUNC(trmv)
154 (const char *uplo, const char *opa, const char *diag, const int *n, const RealScalar *pa, const int *lda,
155  RealScalar *pb, const int *incb) {
156  typedef void (*functype)(int, int, const Scalar *, int, const Scalar *, int, Scalar *, int, const Scalar &);
157  using Eigen::ColMajor;
158  using Eigen::Lower;
159  using Eigen::OnTheLeft;
160  using Eigen::RowMajor;
161  using Eigen::UnitDiag;
162  using Eigen::Upper;
163  static const functype func[16] = {
164  // array index: NOTR | (UP << 2) | (NUNIT << 3)
165  (Eigen::internal::triangular_matrix_vector_product<int, Upper | 0, Scalar, false, Scalar, false, ColMajor>::run),
166  // array index: TR | (UP << 2) | (NUNIT << 3)
167  (Eigen::internal::triangular_matrix_vector_product<int, Lower | 0, Scalar, false, Scalar, false, RowMajor>::run),
168  // array index: ADJ | (UP << 2) | (NUNIT << 3)
169  (Eigen::internal::triangular_matrix_vector_product<int, Lower | 0, Scalar, Conj, Scalar, false, RowMajor>::run),
170  0,
171  // array index: NOTR | (LO << 2) | (NUNIT << 3)
172  (Eigen::internal::triangular_matrix_vector_product<int, Lower | 0, Scalar, false, Scalar, false, ColMajor>::run),
173  // array index: TR | (LO << 2) | (NUNIT << 3)
174  (Eigen::internal::triangular_matrix_vector_product<int, Upper | 0, Scalar, false, Scalar, false, RowMajor>::run),
175  // array index: ADJ | (LO << 2) | (NUNIT << 3)
176  (Eigen::internal::triangular_matrix_vector_product<int, Upper | 0, Scalar, Conj, Scalar, false, RowMajor>::run),
177  0,
178  // array index: NOTR | (UP << 2) | (UNIT << 3)
179  (Eigen::internal::triangular_matrix_vector_product<int, Upper | UnitDiag, Scalar, false, Scalar, false,
180  ColMajor>::run),
181  // array index: TR | (UP << 2) | (UNIT << 3)
182  (Eigen::internal::triangular_matrix_vector_product<int, Lower | UnitDiag, Scalar, false, Scalar, false,
183  RowMajor>::run),
184  // array index: ADJ | (UP << 2) | (UNIT << 3)
185  (Eigen::internal::triangular_matrix_vector_product<int, Lower | UnitDiag, Scalar, Conj, Scalar, false,
186  RowMajor>::run),
187  0,
188  // array index: NOTR | (LO << 2) | (UNIT << 3)
189  (Eigen::internal::triangular_matrix_vector_product<int, Lower | UnitDiag, Scalar, false, Scalar, false,
190  ColMajor>::run),
191  // array index: TR | (LO << 2) | (UNIT << 3)
192  (Eigen::internal::triangular_matrix_vector_product<int, Upper | UnitDiag, Scalar, false, Scalar, false,
193  RowMajor>::run),
194  // array index: ADJ | (LO << 2) | (UNIT << 3)
195  (Eigen::internal::triangular_matrix_vector_product<int, Upper | UnitDiag, Scalar, Conj, Scalar, false,
196  RowMajor>::run),
197  0};
198 
199  const Scalar *a = reinterpret_cast<const Scalar *>(pa);
200  Scalar *b = reinterpret_cast<Scalar *>(pb);
201 
202  int info = 0;
203  if (UPLO(*uplo) == INVALID)
204  info = 1;
205  else if (OP(*opa) == INVALID)
206  info = 2;
207  else if (DIAG(*diag) == INVALID)
208  info = 3;
209  else if (*n < 0)
210  info = 4;
211  else if (*lda < std::max(1, *n))
212  info = 6;
213  else if (*incb == 0)
214  info = 8;
215  if (info) return xerbla_(SCALAR_SUFFIX_UP "TRMV ", &info);
216 
217  if (*n == 0) return;
218 
219  Scalar *actual_b = get_compact_vector(b, *n, *incb);
220  Eigen::Matrix<Scalar, Eigen::Dynamic, 1> res(*n);
221  res.setZero();
222 
223  int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
224  if (code >= 16 || func[code] == 0) return;
225 
226  func[code](*n, *n, a, *lda, actual_b, 1, res.data(), 1, Scalar(1));
227 
228  copy_back(res.data(), b, *n, *incb);
229  if (actual_b != b) delete[] actual_b;
230 }
231 
239 EIGEN_BLAS_FUNC(gbmv)
240 (char *trans, int *m, int *n, int *kl, int *ku, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *px, int *incx,
241  RealScalar *pbeta, RealScalar *py, int *incy) {
242  const Scalar *a = reinterpret_cast<const Scalar *>(pa);
243  const Scalar *x = reinterpret_cast<const Scalar *>(px);
244  Scalar *y = reinterpret_cast<Scalar *>(py);
245  Scalar alpha = *reinterpret_cast<const Scalar *>(palpha);
246  Scalar beta = *reinterpret_cast<const Scalar *>(pbeta);
247  int coeff_rows = *kl + *ku + 1;
248 
249  int info = 0;
250  if (OP(*trans) == INVALID)
251  info = 1;
252  else if (*m < 0)
253  info = 2;
254  else if (*n < 0)
255  info = 3;
256  else if (*kl < 0)
257  info = 4;
258  else if (*ku < 0)
259  info = 5;
260  else if (*lda < coeff_rows)
261  info = 8;
262  else if (*incx == 0)
263  info = 10;
264  else if (*incy == 0)
265  info = 13;
266  if (info) return xerbla_(SCALAR_SUFFIX_UP "GBMV ", &info);
267 
268  if (*m == 0 || *n == 0 || (alpha == Scalar(0) && beta == Scalar(1))) return;
269 
270  int actual_m = *m;
271  int actual_n = *n;
272  if (OP(*trans) != NOTR) std::swap(actual_m, actual_n);
273 
274  const Scalar *actual_x = get_compact_vector(x, actual_n, *incx);
275  Scalar *actual_y = get_compact_vector(y, actual_m, *incy);
276 
277  if (beta != Scalar(1)) {
278  if (beta == Scalar(0))
279  make_vector(actual_y, actual_m).setZero();
280  else
281  make_vector(actual_y, actual_m) *= beta;
282  }
283 
284  ConstMatrixType mat_coeffs(a, coeff_rows, *n, *lda);
285 
286  int nb = std::min(*n, (*m) + (*ku));
287  for (int j = 0; j < nb; ++j) {
288  int start = std::max(0, j - *ku);
289  int end = std::min((*m) - 1, j + *kl);
290  int len = end - start + 1;
291  int offset = (*ku) - j + start;
292  if (OP(*trans) == NOTR)
293  make_vector(actual_y + start, len) += (alpha * actual_x[j]) * mat_coeffs.col(j).segment(offset, len);
294  else if (OP(*trans) == TR)
295  actual_y[j] +=
296  alpha * (mat_coeffs.col(j).segment(offset, len).transpose() * make_vector(actual_x + start, len)).value();
297  else
298  actual_y[j] +=
299  alpha * (mat_coeffs.col(j).segment(offset, len).adjoint() * make_vector(actual_x + start, len)).value();
300  }
301 
302  if (actual_x != x) delete[] actual_x;
303  if (actual_y != y) delete[] copy_back(actual_y, y, actual_m, *incy);
304 }
305 
306 #if 0
314 EIGEN_BLAS_FUNC(tbmv)(char *uplo, char *opa, char *diag, int *n, int *k, RealScalar *pa, int *lda, RealScalar *px, int *incx)
315 {
316  Scalar* a = reinterpret_cast<Scalar*>(pa);
317  Scalar* x = reinterpret_cast<Scalar*>(px);
318  int coeff_rows = *k + 1;
319 
320  int info = 0;
321  if(UPLO(*uplo)==INVALID) info = 1;
322  else if(OP(*opa)==INVALID) info = 2;
323  else if(DIAG(*diag)==INVALID) info = 3;
324  else if(*n<0) info = 4;
325  else if(*k<0) info = 5;
326  else if(*lda<coeff_rows) info = 7;
327  else if(*incx==0) info = 9;
328  if(info)
329  return xerbla_(SCALAR_SUFFIX_UP"TBMV ",&info,6);
330 
331  if(*n==0) return;
332 
333  int actual_n = *n;
334 
335  Scalar* actual_x = get_compact_vector(x,actual_n,*incx);
336 
337  MatrixType mat_coeffs(a,coeff_rows,*n,*lda);
338 
339  int ku = UPLO(*uplo)==UPPER ? *k : 0;
340  int kl = UPLO(*uplo)==LOWER ? *k : 0;
341 
342  for(int j=0; j<*n; ++j)
343  {
344  int start = std::max(0,j - ku);
345  int end = std::min((*m)-1,j + kl);
346  int len = end - start + 1;
347  int offset = (ku) - j + start;
348 
349  if(OP(*trans)==NOTR)
350  make_vector(actual_y+start,len) += (alpha*actual_x[j]) * mat_coeffs.col(j).segment(offset,len);
351  else if(OP(*trans)==TR)
352  actual_y[j] += alpha * ( mat_coeffs.col(j).segment(offset,len).transpose() * make_vector(actual_x+start,len) ).value();
353  else
354  actual_y[j] += alpha * ( mat_coeffs.col(j).segment(offset,len).adjoint() * make_vector(actual_x+start,len) ).value();
355  }
356 
357  if(actual_x!=x) delete[] actual_x;
358  if(actual_y!=y) delete[] copy_back(actual_y,y,actual_m,*incy);
359 }
360 #endif
361 
373 EIGEN_BLAS_FUNC(tbsv)
374 (char *uplo, char *op, char *diag, int *n, int *k, RealScalar *pa, int *lda, RealScalar *px, int *incx) {
375  typedef void (*functype)(int, int, const Scalar *, int, Scalar *);
376  using Eigen::ColMajor;
377  using Eigen::Lower;
378  using Eigen::OnTheLeft;
379  using Eigen::RowMajor;
380  using Eigen::UnitDiag;
381  using Eigen::Upper;
382  static const functype func[16] = {
383  // array index: NOTR | (UP << 2) | (NUNIT << 3)
384  (Eigen::internal::band_solve_triangular_selector<int, Upper | 0, Scalar, false, Scalar, ColMajor>::run),
385  // array index: TR | (UP << 2) | (NUNIT << 3)
386  (Eigen::internal::band_solve_triangular_selector<int, Lower | 0, Scalar, false, Scalar, RowMajor>::run),
387  // array index: ADJ | (UP << 2) | (NUNIT << 3)
388  (Eigen::internal::band_solve_triangular_selector<int, Lower | 0, Scalar, Conj, Scalar, RowMajor>::run),
389  0,
390  // array index: NOTR | (LO << 2) | (NUNIT << 3)
391  (Eigen::internal::band_solve_triangular_selector<int, Lower | 0, Scalar, false, Scalar, ColMajor>::run),
392  // array index: TR | (LO << 2) | (NUNIT << 3)
393  (Eigen::internal::band_solve_triangular_selector<int, Upper | 0, Scalar, false, Scalar, RowMajor>::run),
394  // array index: ADJ | (LO << 2) | (NUNIT << 3)
395  (Eigen::internal::band_solve_triangular_selector<int, Upper | 0, Scalar, Conj, Scalar, RowMajor>::run),
396  0,
397  // array index: NOTR | (UP << 2) | (UNIT << 3)
398  (Eigen::internal::band_solve_triangular_selector<int, Upper | UnitDiag, Scalar, false, Scalar, ColMajor>::run),
399  // array index: TR | (UP << 2) | (UNIT << 3)
400  (Eigen::internal::band_solve_triangular_selector<int, Lower | UnitDiag, Scalar, false, Scalar, RowMajor>::run),
401  // array index: ADJ | (UP << 2) | (UNIT << 3)
402  (Eigen::internal::band_solve_triangular_selector<int, Lower | UnitDiag, Scalar, Conj, Scalar, RowMajor>::run),
403  0,
404  // array index: NOTR | (LO << 2) | (UNIT << 3)
405  (Eigen::internal::band_solve_triangular_selector<int, Lower | UnitDiag, Scalar, false, Scalar, ColMajor>::run),
406  // array index: TR | (LO << 2) | (UNIT << 3)
407  (Eigen::internal::band_solve_triangular_selector<int, Upper | UnitDiag, Scalar, false, Scalar, RowMajor>::run),
408  // array index: ADJ | (LO << 2) | (UNIT << 3)
409  (Eigen::internal::band_solve_triangular_selector<int, Upper | UnitDiag, Scalar, Conj, Scalar, RowMajor>::run),
410  0,
411  };
412 
413  Scalar *a = reinterpret_cast<Scalar *>(pa);
414  Scalar *x = reinterpret_cast<Scalar *>(px);
415  int coeff_rows = *k + 1;
416 
417  int info = 0;
418  if (UPLO(*uplo) == INVALID)
419  info = 1;
420  else if (OP(*op) == INVALID)
421  info = 2;
422  else if (DIAG(*diag) == INVALID)
423  info = 3;
424  else if (*n < 0)
425  info = 4;
426  else if (*k < 0)
427  info = 5;
428  else if (*lda < coeff_rows)
429  info = 7;
430  else if (*incx == 0)
431  info = 9;
432  if (info) return xerbla_(SCALAR_SUFFIX_UP "TBSV ", &info);
433 
434  if (*n == 0 || (*k == 0 && DIAG(*diag) == UNIT)) return;
435 
436  int actual_n = *n;
437 
438  Scalar *actual_x = get_compact_vector(x, actual_n, *incx);
439 
440  int code = OP(*op) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
441  if (code >= 16 || func[code] == 0) return;
442 
443  func[code](*n, *k, a, *lda, actual_x);
444 
445  if (actual_x != x) delete[] copy_back(actual_x, x, actual_n, *incx);
446 }
447 
455 EIGEN_BLAS_FUNC(tpmv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pap, RealScalar *px, int *incx) {
456  typedef void (*functype)(int, const Scalar *, const Scalar *, Scalar *, Scalar);
457  using Eigen::ColMajor;
458  using Eigen::Lower;
459  using Eigen::OnTheLeft;
460  using Eigen::RowMajor;
461  using Eigen::UnitDiag;
462  using Eigen::Upper;
463  static const functype func[16] = {
464  // array index: NOTR | (UP << 2) | (NUNIT << 3)
465  (Eigen::internal::packed_triangular_matrix_vector_product<int, Upper | 0, Scalar, false, Scalar, false,
466  ColMajor>::run),
467  // array index: TR | (UP << 2) | (NUNIT << 3)
468  (Eigen::internal::packed_triangular_matrix_vector_product<int, Lower | 0, Scalar, false, Scalar, false,
469  RowMajor>::run),
470  // array index: ADJ | (UP << 2) | (NUNIT << 3)
471  (Eigen::internal::packed_triangular_matrix_vector_product<int, Lower | 0, Scalar, Conj, Scalar, false,
472  RowMajor>::run),
473  0,
474  // array index: NOTR | (LO << 2) | (NUNIT << 3)
475  (Eigen::internal::packed_triangular_matrix_vector_product<int, Lower | 0, Scalar, false, Scalar, false,
476  ColMajor>::run),
477  // array index: TR | (LO << 2) | (NUNIT << 3)
478  (Eigen::internal::packed_triangular_matrix_vector_product<int, Upper | 0, Scalar, false, Scalar, false,
479  RowMajor>::run),
480  // array index: ADJ | (LO << 2) | (NUNIT << 3)
481  (Eigen::internal::packed_triangular_matrix_vector_product<int, Upper | 0, Scalar, Conj, Scalar, false,
482  RowMajor>::run),
483  0,
484  // array index: NOTR | (UP << 2) | (UNIT << 3)
485  (Eigen::internal::packed_triangular_matrix_vector_product<int, Upper | UnitDiag, Scalar, false, Scalar, false,
486  ColMajor>::run),
487  // array index: TR | (UP << 2) | (UNIT << 3)
488  (Eigen::internal::packed_triangular_matrix_vector_product<int, Lower | UnitDiag, Scalar, false, Scalar, false,
489  RowMajor>::run),
490  // array index: ADJ | (UP << 2) | (UNIT << 3)
491  (Eigen::internal::packed_triangular_matrix_vector_product<int, Lower | UnitDiag, Scalar, Conj, Scalar, false,
492  RowMajor>::run),
493  0,
494  // array index: NOTR | (LO << 2) | (UNIT << 3)
495  (Eigen::internal::packed_triangular_matrix_vector_product<int, Lower | UnitDiag, Scalar, false, Scalar, false,
496  ColMajor>::run),
497  // array index: TR | (LO << 2) | (UNIT << 3)
498  (Eigen::internal::packed_triangular_matrix_vector_product<int, Upper | UnitDiag, Scalar, false, Scalar, false,
499  RowMajor>::run),
500  // array index: ADJ | (LO << 2) | (UNIT << 3)
501  (Eigen::internal::packed_triangular_matrix_vector_product<int, Upper | UnitDiag, Scalar, Conj, Scalar, false,
502  RowMajor>::run),
503  0};
504 
505  Scalar *ap = reinterpret_cast<Scalar *>(pap);
506  Scalar *x = reinterpret_cast<Scalar *>(px);
507 
508  int info = 0;
509  if (UPLO(*uplo) == INVALID)
510  info = 1;
511  else if (OP(*opa) == INVALID)
512  info = 2;
513  else if (DIAG(*diag) == INVALID)
514  info = 3;
515  else if (*n < 0)
516  info = 4;
517  else if (*incx == 0)
518  info = 7;
519  if (info) return xerbla_(SCALAR_SUFFIX_UP "TPMV ", &info);
520 
521  if (*n == 0) return;
522 
523  Scalar *actual_x = get_compact_vector(x, *n, *incx);
524  Eigen::Matrix<Scalar, Eigen::Dynamic, 1> res(*n);
525  res.setZero();
526 
527  int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
528  if (code >= 16 || func[code] == 0) return;
529 
530  func[code](*n, ap, actual_x, res.data(), Scalar(1));
531 
532  copy_back(res.data(), x, *n, *incx);
533  if (actual_x != x) delete[] actual_x;
534 }
535 
546 EIGEN_BLAS_FUNC(tpsv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pap, RealScalar *px, int *incx) {
547  typedef void (*functype)(int, const Scalar *, Scalar *);
548  using Eigen::ColMajor;
549  using Eigen::Lower;
550  using Eigen::OnTheLeft;
551  using Eigen::RowMajor;
552  using Eigen::UnitDiag;
553  using Eigen::Upper;
554  static const functype func[16] = {
555  // array index: NOTR | (UP << 2) | (NUNIT << 3)
556  (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Upper | 0, false,
557  ColMajor>::run),
558  // array index: TR | (UP << 2) | (NUNIT << 3)
559  (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Lower | 0, false,
560  RowMajor>::run),
561  // array index: ADJ | (UP << 2) | (NUNIT << 3)
562  (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Lower | 0, Conj, RowMajor>::run),
563  0,
564  // array index: NOTR | (LO << 2) | (NUNIT << 3)
565  (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Lower | 0, false,
566  ColMajor>::run),
567  // array index: TR | (LO << 2) | (NUNIT << 3)
568  (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Upper | 0, false,
569  RowMajor>::run),
570  // array index: ADJ | (LO << 2) | (NUNIT << 3)
571  (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Upper | 0, Conj, RowMajor>::run),
572  0,
573  // array index: NOTR | (UP << 2) | (UNIT << 3)
574  (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Upper | UnitDiag, false,
575  ColMajor>::run),
576  // array index: TR | (UP << 2) | (UNIT << 3)
577  (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Lower | UnitDiag, false,
578  RowMajor>::run),
579  // array index: ADJ | (UP << 2) | (UNIT << 3)
580  (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Lower | UnitDiag, Conj,
581  RowMajor>::run),
582  0,
583  // array index: NOTR | (LO << 2) | (UNIT << 3)
584  (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Lower | UnitDiag, false,
585  ColMajor>::run),
586  // array index: TR | (LO << 2) | (UNIT << 3)
587  (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Upper | UnitDiag, false,
588  RowMajor>::run),
589  // array index: ADJ | (LO << 2) | (UNIT << 3)
590  (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Upper | UnitDiag, Conj,
591  RowMajor>::run),
592  0};
593 
594  Scalar *ap = reinterpret_cast<Scalar *>(pap);
595  Scalar *x = reinterpret_cast<Scalar *>(px);
596 
597  int info = 0;
598  if (UPLO(*uplo) == INVALID)
599  info = 1;
600  else if (OP(*opa) == INVALID)
601  info = 2;
602  else if (DIAG(*diag) == INVALID)
603  info = 3;
604  else if (*n < 0)
605  info = 4;
606  else if (*incx == 0)
607  info = 7;
608  if (info) return xerbla_(SCALAR_SUFFIX_UP "TPSV ", &info);
609 
610  Scalar *actual_x = get_compact_vector(x, *n, *incx);
611 
612  int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
613  func[code](*n, ap, actual_x);
614 
615  if (actual_x != x) delete[] copy_back(actual_x, x, *n, *incx);
616 }
res
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
Definition: PartialRedux_count.cpp:3
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:45
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:46
MatrixType
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
SCALAR_SUFFIX_UP
#define SCALAR_SUFFIX_UP
Definition: blas/complex_double.cpp:12
Eigen::Map
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:96
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: Eigen/Eigen/src/Core/Matrix.h:186
Eigen::internal::const_blas_data_mapper
Definition: BlasUtil.h:443
UNIT
#define UNIT
Definition: common.h:50
make_vector
Eigen::Map< Eigen::Matrix< T, Eigen::Dynamic, 1 >, 0, Eigen::InnerStride< Eigen::Dynamic > > make_vector(T *data, int size, int incr)
Definition: common.h:99
Conj
@ Conj
Definition: common.h:73
get_compact_vector
T * get_compact_vector(T *x, int n, int incx)
Definition: common.h:124
TR
#define TR
Definition: common.h:40
OP
#define OP(X)
Definition: common.h:54
NOTR
#define NOTR
Definition: common.h:39
INVALID
#define INVALID
Definition: common.h:52
UPLO
#define UPLO(X)
Definition: common.h:59
DIAG
#define DIAG(X)
Definition: common.h:61
min
#define min(a, b)
Definition: datatypes.h:22
max
#define max(a, b)
Definition: datatypes.h:23
gemv
EIGEN_DONT_INLINE void gemv(const Mat &A, const Vec &B, Vec &C)
Definition: gemv.cpp:3
Eigen::placeholders::end
static constexpr lastp1_t end
Definition: IndexedViewHelper.h:79
Eigen::UnitDiag
@ UnitDiag
Definition: Constants.h:215
Eigen::Lower
@ Lower
Definition: Constants.h:211
Eigen::Upper
@ Upper
Definition: Constants.h:213
Eigen::ColMajor
@ ColMajor
Definition: Constants.h:318
Eigen::RowMajor
@ RowMajor
Definition: Constants.h:320
Eigen::OnTheLeft
@ OnTheLeft
Definition: Constants.h:331
int
return int(ret)+1
EIGEN_BLAS_FUNC
EIGEN_BLAS_FUNC(EIGEN_CAT(REAL_SCALAR_SUFFIX, scal))(int *n
swap
EIGEN_BLAS_FUNC() swap(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
Definition: level1_impl.h:117
ap
Scalar * ap
Definition: level2_cplx_impl.h:161
pap
char int RealScalar RealScalar int RealScalar int RealScalar * pap
Definition: level2_cplx_impl.h:150
pbeta
const char const int const int const RealScalar const RealScalar const int const RealScalar const int const RealScalar * pbeta
Definition: level2_impl.h:28
actual_m
actual_m
Definition: level2_impl.h:82
incy
char int int int int RealScalar RealScalar int RealScalar int RealScalar RealScalar int * incy
Definition: level2_impl.h:241
px
char int int int int RealScalar RealScalar int RealScalar * px
Definition: level2_impl.h:240
palpha
const char const int const int const RealScalar * palpha
Definition: level2_impl.h:27
coeff_rows
int coeff_rows
Definition: level2_impl.h:247
y
Scalar * y
Definition: level2_impl.h:244
ku
char int int int int * ku
Definition: level2_impl.h:240
code
int code
Definition: level2_impl.h:147
actual_b
if(code !=NOTR) std const Scalar * actual_b
Definition: level2_impl.h:67
actual_n
actual_n
Definition: level2_impl.h:445
info
int info
Definition: level2_impl.h:45
a
const Scalar * a
Definition: level2_impl.h:38
trans
char * trans
Definition: level2_impl.h:240
b
const Scalar * b
Definition: level2_impl.h:39
n
const char const int const int * n
Definition: level2_impl.h:27
copy_back
copy_back(res.data(), b, *n, *incb)
actual_c
Scalar * actual_c
Definition: level2_impl.h:68
uplo
const char * uplo
Definition: level2_impl.h:86
x
const Scalar * x
Definition: level2_impl.h:243
lda
const char const int const int const RealScalar const RealScalar const int * lda
Definition: level2_impl.h:27
incb
const char const int const int const RealScalar const RealScalar const int const RealScalar const int * incb
Definition: level2_impl.h:28
diag
const char const char const char * diag
Definition: level2_impl.h:86
py
char int int int int RealScalar RealScalar int RealScalar int RealScalar RealScalar * py
Definition: level2_impl.h:241
actual_y
Scalar * actual_y
Definition: level2_impl.h:275
incc
const char const int const int const RealScalar const RealScalar const int const RealScalar const int const RealScalar RealScalar const int * incc
Definition: level2_impl.h:28
kl
char int int int * kl
Definition: level2_impl.h:240
alpha
Scalar alpha
Definition: level2_impl.h:41
beta
Scalar beta
Definition: level2_impl.h:42
mat_coeffs
ConstMatrixType mat_coeffs(a, coeff_rows, *n, *lda)
m
const char const int * m
Definition: level2_impl.h:27
nb
int nb
Definition: level2_impl.h:286
pb
const char const int const int const RealScalar const RealScalar const int const RealScalar * pb
Definition: level2_impl.h:28
opa
const char * opa
Definition: level2_impl.h:27
actual_x
const Scalar * actual_x
Definition: level2_impl.h:274
k
char char char int int * k
Definition: level2_impl.h:374
op
char char * op
Definition: level2_impl.h:374
pc
const char const int const int const RealScalar const RealScalar const int const RealScalar const int const RealScalar RealScalar * pc
Definition: level2_impl.h:28
incx
char int int int int RealScalar RealScalar int RealScalar int * incx
Definition: level2_impl.h:240
c
Scalar * c
Definition: level2_impl.h:40
tpsv
EIGEN_BLAS_FUNC() tpsv(char *uplo, char *opa, char *diag, int *n, RealScalar *pap, RealScalar *px, int *incx)
Definition: level2_impl.h:546
pa
const char const int const int const RealScalar const RealScalar * pa
Definition: level2_impl.h:27
tpmv
EIGEN_BLAS_FUNC() tpmv(char *uplo, char *opa, char *diag, int *n, RealScalar *pap, RealScalar *px, int *incx)
Definition: level2_impl.h:455
functype
void(* functype)(DenseIndex, DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, DenseIndex, Scalar, Eigen::internal::level3_blocking< Scalar, Scalar > &, Eigen::internal::GemmParallelInfo< DenseIndex > *)
Definition: level3_impl.h:22
Eigen::value
squared absolute value
Definition: GlobalFunctions.h:87
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:83
oomph::CumulativeTimings::start
void start(const unsigned &i)
(Re-)start i-th timer
Definition: oomph_utilities.cc:243
Eigen::internal::band_solve_triangular_selector
Definition: BandTriangularSolver.h:20
Eigen::internal::general_matrix_vector_product
Definition: BlasUtil.h:42
Eigen::internal::packed_triangular_matrix_vector_product
Definition: PackedTriangularMatrixVector.h:18
Eigen::internal::packed_triangular_solve_vector
Definition: PackedTriangularSolverVector.h:17
Eigen::internal::triangular_matrix_vector_product
Definition: TriangularMatrixVector.h:22
Eigen::internal::triangular_solve_vector
Definition: SolveTriangular.h:23
func
Definition: benchGeometry.cpp:21
general_matrix_vector_product_wrapper
Definition: level2_impl.h:13
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:14
trmv
EIGEN_DONT_INLINE void trmv(const Mat &A, const Vec &B, Vec &C)
Definition: trmv_lo.cpp:3
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2
run
void run(const string &dir_name, LinearSolver *linear_solver_pt, const unsigned nel_1d, bool mess_up_order)
Definition: two_d_poisson_compare_solvers.cc:317
xerbla_
EIGEN_WEAK_LINKING void xerbla_(const char *msg, int *info)
Definition: xerbla.cpp:14