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 {
rungeneral_matrix_vector_product_wrapper15 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
EIGEN_BLAS_FUNC(gemv)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
EIGEN_BLAS_FUNC(trsv)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
EIGEN_BLAS_FUNC(trmv)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
211 /** GBMV performs one of the matrix-vector operations
212 *
213 * y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y,
214 *
215 * where alpha and beta are scalars, x and y are vectors and A is an
216 * m by n band matrix, with kl sub-diagonals and ku super-diagonals.
217 */
EIGEN_BLAS_FUNC(gbmv)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 /** TBMV performs one of the matrix-vector operations
282 *
283 * x := A*x, or x := A'*x,
284 *
285 * where x is an n element vector and A is an n by n unit, or non-unit,
286 * upper or lower triangular band matrix, with ( k + 1 ) diagonals.
287 */
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
339 /** DTBSV solves one of the systems of equations
340 *
341 * A*x = b, or A'*x = b,
342 *
343 * where b and x are n element vectors and A is an n by n unit, or
344 * non-unit, upper or lower triangular band matrix, with ( k + 1 )
345 * diagonals.
346 *
347 * No test for singularity or near-singularity is included in this
348 * routine. Such tests must be performed before calling this routine.
349 */
EIGEN_BLAS_FUNC(tbsv)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
417 /** DTPMV performs one of the matrix-vector operations
418 *
419 * x := A*x, or x := A'*x,
420 *
421 * where x is an n element vector and A is an n by n unit, or non-unit,
422 * upper or lower triangular matrix, supplied in packed form.
423 */
EIGEN_BLAS_FUNC(tpmv)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
489 /** DTPSV solves one of the systems of equations
490 *
491 * A*x = b, or A'*x = b,
492 *
493 * where b and x are n element vectors and A is an n by n unit, or
494 * non-unit, upper or lower triangular matrix, supplied in packed form.
495 *
496 * No test for singularity or near-singularity is included in this
497 * routine. Such tests must be performed before calling this routine.
498 */
EIGEN_BLAS_FUNC(tpsv)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 }
554