RavEngine
Loading...
Searching...
No Matches
level2_impl.h
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
12template<typename Index, typename Scalar, int StorageOrder, bool ConjugateLhs, bool ConjugateRhs>
13struct 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
26int 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
86int 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
145int 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
218int 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
288int 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
350int 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
424int 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
499int 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}
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 Matrix.h:180
Eigen::internal::const_blas_data_mapper
Definition BlasUtil.h:389
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition Meta.h:74
Eigen::internal::general_matrix_vector_product
Definition BlasUtil.h:40
general_matrix_vector_product_wrapper
Definition level2_impl.h:14
internal::band_solve_triangular_selector
Definition BandTriangularSolver.h:19
internal::packed_triangular_matrix_vector_product
Definition PackedTriangularMatrixVector.h:16
internal::packed_triangular_solve_vector
Definition PackedTriangularSolverVector.h:16