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
EIGEN_BLAS_FUNC(gemv)12 int EIGEN_BLAS_FUNC(gemv)(char *opa, int *m, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *incb, RealScalar *pbeta, RealScalar *pc, int *incc)
13 {
14 typedef void (*functype)(int, int, const Scalar *, int, const Scalar *, int , Scalar *, int, Scalar);
15 static functype func[4];
16
17 static bool init = false;
18 if(!init)
19 {
20 for(int k=0; k<4; ++k)
21 func[k] = 0;
22
23 func[NOTR] = (internal::general_matrix_vector_product<int,Scalar,ColMajor,false,Scalar,false>::run);
24 func[TR ] = (internal::general_matrix_vector_product<int,Scalar,RowMajor,false,Scalar,false>::run);
25 func[ADJ ] = (internal::general_matrix_vector_product<int,Scalar,RowMajor,Conj, Scalar,false>::run);
26
27 init = true;
28 }
29
30 Scalar* a = reinterpret_cast<Scalar*>(pa);
31 Scalar* b = reinterpret_cast<Scalar*>(pb);
32 Scalar* c = reinterpret_cast<Scalar*>(pc);
33 Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
34 Scalar beta = *reinterpret_cast<Scalar*>(pbeta);
35
36 // check arguments
37 int info = 0;
38 if(OP(*opa)==INVALID) info = 1;
39 else if(*m<0) info = 2;
40 else if(*n<0) info = 3;
41 else if(*lda<std::max(1,*m)) info = 6;
42 else if(*incb==0) info = 8;
43 else if(*incc==0) info = 11;
44 if(info)
45 return xerbla_(SCALAR_SUFFIX_UP"GEMV ",&info,6);
46
47 if(*m==0 || *n==0 || (alpha==Scalar(0) && beta==Scalar(1)))
48 return 0;
49
50 int actual_m = *m;
51 int actual_n = *n;
52 int code = OP(*opa);
53 if(code!=NOTR)
54 std::swap(actual_m,actual_n);
55
56 Scalar* actual_b = get_compact_vector(b,actual_n,*incb);
57 Scalar* actual_c = get_compact_vector(c,actual_m,*incc);
58
59 if(beta!=Scalar(1))
60 {
61 if(beta==Scalar(0)) vector(actual_c, actual_m).setZero();
62 else vector(actual_c, actual_m) *= beta;
63 }
64
65 if(code>=4 || func[code]==0)
66 return 0;
67
68 func[code](actual_m, actual_n, a, *lda, actual_b, 1, actual_c, 1, alpha);
69
70 if(actual_b!=b) delete[] actual_b;
71 if(actual_c!=c) delete[] copy_back(actual_c,c,actual_m,*incc);
72
73 return 1;
74 }
75
EIGEN_BLAS_FUNC(trsv)76 int EIGEN_BLAS_FUNC(trsv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pa, int *lda, RealScalar *pb, int *incb)
77 {
78 typedef void (*functype)(int, const Scalar *, int, Scalar *);
79 static functype func[16];
80
81 static bool init = false;
82 if(!init)
83 {
84 for(int k=0; k<16; ++k)
85 func[k] = 0;
86
87 func[NOTR | (UP << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0, false,ColMajor>::run);
88 func[TR | (UP << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0, false,RowMajor>::run);
89 func[ADJ | (UP << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0, Conj, RowMajor>::run);
90
91 func[NOTR | (LO << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0, false,ColMajor>::run);
92 func[TR | (LO << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0, false,RowMajor>::run);
93 func[ADJ | (LO << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0, Conj, RowMajor>::run);
94
95 func[NOTR | (UP << 2) | (UNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,ColMajor>::run);
96 func[TR | (UP << 2) | (UNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,RowMajor>::run);
97 func[ADJ | (UP << 2) | (UNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,Conj, RowMajor>::run);
98
99 func[NOTR | (LO << 2) | (UNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,ColMajor>::run);
100 func[TR | (LO << 2) | (UNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,RowMajor>::run);
101 func[ADJ | (LO << 2) | (UNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,Conj, RowMajor>::run);
102
103 init = true;
104 }
105
106 Scalar* a = reinterpret_cast<Scalar*>(pa);
107 Scalar* b = reinterpret_cast<Scalar*>(pb);
108
109 int info = 0;
110 if(UPLO(*uplo)==INVALID) info = 1;
111 else if(OP(*opa)==INVALID) info = 2;
112 else if(DIAG(*diag)==INVALID) info = 3;
113 else if(*n<0) info = 4;
114 else if(*lda<std::max(1,*n)) info = 6;
115 else if(*incb==0) info = 8;
116 if(info)
117 return xerbla_(SCALAR_SUFFIX_UP"TRSV ",&info,6);
118
119 Scalar* actual_b = get_compact_vector(b,*n,*incb);
120
121 int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
122 func[code](*n, a, *lda, actual_b);
123
124 if(actual_b!=b) delete[] copy_back(actual_b,b,*n,*incb);
125
126 return 0;
127 }
128
129
130
EIGEN_BLAS_FUNC(trmv)131 int EIGEN_BLAS_FUNC(trmv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pa, int *lda, RealScalar *pb, int *incb)
132 {
133 typedef void (*functype)(int, int, const Scalar *, int, const Scalar *, int, Scalar *, int, const Scalar&);
134 static functype func[16];
135
136 static bool init = false;
137 if(!init)
138 {
139 for(int k=0; k<16; ++k)
140 func[k] = 0;
141
142 func[NOTR | (UP << 2) | (NUNIT << 3)] = (internal::triangular_matrix_vector_product<int,Upper|0, Scalar,false,Scalar,false,ColMajor>::run);
143 func[TR | (UP << 2) | (NUNIT << 3)] = (internal::triangular_matrix_vector_product<int,Lower|0, Scalar,false,Scalar,false,RowMajor>::run);
144 func[ADJ | (UP << 2) | (NUNIT << 3)] = (internal::triangular_matrix_vector_product<int,Lower|0, Scalar,Conj, Scalar,false,RowMajor>::run);
145
146 func[NOTR | (LO << 2) | (NUNIT << 3)] = (internal::triangular_matrix_vector_product<int,Lower|0, Scalar,false,Scalar,false,ColMajor>::run);
147 func[TR | (LO << 2) | (NUNIT << 3)] = (internal::triangular_matrix_vector_product<int,Upper|0, Scalar,false,Scalar,false,RowMajor>::run);
148 func[ADJ | (LO << 2) | (NUNIT << 3)] = (internal::triangular_matrix_vector_product<int,Upper|0, Scalar,Conj, Scalar,false,RowMajor>::run);
149
150 func[NOTR | (UP << 2) | (UNIT << 3)] = (internal::triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run);
151 func[TR | (UP << 2) | (UNIT << 3)] = (internal::triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run);
152 func[ADJ | (UP << 2) | (UNIT << 3)] = (internal::triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run);
153
154 func[NOTR | (LO << 2) | (UNIT << 3)] = (internal::triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run);
155 func[TR | (LO << 2) | (UNIT << 3)] = (internal::triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run);
156 func[ADJ | (LO << 2) | (UNIT << 3)] = (internal::triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run);
157
158 init = true;
159 }
160
161 Scalar* a = reinterpret_cast<Scalar*>(pa);
162 Scalar* b = reinterpret_cast<Scalar*>(pb);
163
164 int info = 0;
165 if(UPLO(*uplo)==INVALID) info = 1;
166 else if(OP(*opa)==INVALID) info = 2;
167 else if(DIAG(*diag)==INVALID) info = 3;
168 else if(*n<0) info = 4;
169 else if(*lda<std::max(1,*n)) info = 6;
170 else if(*incb==0) info = 8;
171 if(info)
172 return xerbla_(SCALAR_SUFFIX_UP"TRMV ",&info,6);
173
174 if(*n==0)
175 return 1;
176
177 Scalar* actual_b = get_compact_vector(b,*n,*incb);
178 Matrix<Scalar,Dynamic,1> res(*n);
179 res.setZero();
180
181 int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
182 if(code>=16 || func[code]==0)
183 return 0;
184
185 func[code](*n, *n, a, *lda, actual_b, 1, res.data(), 1, Scalar(1));
186
187 copy_back(res.data(),b,*n,*incb);
188 if(actual_b!=b) delete[] actual_b;
189
190 return 1;
191 }
192
193 /** GBMV performs one of the matrix-vector operations
194 *
195 * y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y,
196 *
197 * where alpha and beta are scalars, x and y are vectors and A is an
198 * m by n band matrix, with kl sub-diagonals and ku super-diagonals.
199 */
EIGEN_BLAS_FUNC(gbmv)200 int EIGEN_BLAS_FUNC(gbmv)(char *trans, int *m, int *n, int *kl, int *ku, RealScalar *palpha, RealScalar *pa, int *lda,
201 RealScalar *px, int *incx, RealScalar *pbeta, RealScalar *py, int *incy)
202 {
203 Scalar* a = reinterpret_cast<Scalar*>(pa);
204 Scalar* x = reinterpret_cast<Scalar*>(px);
205 Scalar* y = reinterpret_cast<Scalar*>(py);
206 Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
207 Scalar beta = *reinterpret_cast<Scalar*>(pbeta);
208 int coeff_rows = *kl+*ku+1;
209
210 int info = 0;
211 if(OP(*trans)==INVALID) info = 1;
212 else if(*m<0) info = 2;
213 else if(*n<0) info = 3;
214 else if(*kl<0) info = 4;
215 else if(*ku<0) info = 5;
216 else if(*lda<coeff_rows) info = 8;
217 else if(*incx==0) info = 10;
218 else if(*incy==0) info = 13;
219 if(info)
220 return xerbla_(SCALAR_SUFFIX_UP"GBMV ",&info,6);
221
222 if(*m==0 || *n==0 || (alpha==Scalar(0) && beta==Scalar(1)))
223 return 0;
224
225 int actual_m = *m;
226 int actual_n = *n;
227 if(OP(*trans)!=NOTR)
228 std::swap(actual_m,actual_n);
229
230 Scalar* actual_x = get_compact_vector(x,actual_n,*incx);
231 Scalar* actual_y = get_compact_vector(y,actual_m,*incy);
232
233 if(beta!=Scalar(1))
234 {
235 if(beta==Scalar(0)) vector(actual_y, actual_m).setZero();
236 else vector(actual_y, actual_m) *= beta;
237 }
238
239 MatrixType mat_coeffs(a,coeff_rows,*n,*lda);
240
241 int nb = std::min(*n,(*m)+(*ku));
242 for(int j=0; j<nb; ++j)
243 {
244 int start = std::max(0,j - *ku);
245 int end = std::min((*m)-1,j + *kl);
246 int len = end - start + 1;
247 int offset = (*ku) - j + start;
248 if(OP(*trans)==NOTR)
249 vector(actual_y+start,len) += (alpha*actual_x[j]) * mat_coeffs.col(j).segment(offset,len);
250 else if(OP(*trans)==TR)
251 actual_y[j] += alpha * ( mat_coeffs.col(j).segment(offset,len).transpose() * vector(actual_x+start,len) ).value();
252 else
253 actual_y[j] += alpha * ( mat_coeffs.col(j).segment(offset,len).adjoint() * vector(actual_x+start,len) ).value();
254 }
255
256 if(actual_x!=x) delete[] actual_x;
257 if(actual_y!=y) delete[] copy_back(actual_y,y,actual_m,*incy);
258
259 return 0;
260 }
261
262 #if 0
263 /** TBMV performs one of the matrix-vector operations
264 *
265 * x := A*x, or x := A'*x,
266 *
267 * where x is an n element vector and A is an n by n unit, or non-unit,
268 * upper or lower triangular band matrix, with ( k + 1 ) diagonals.
269 */
270 int EIGEN_BLAS_FUNC(tbmv)(char *uplo, char *opa, char *diag, int *n, int *k, RealScalar *pa, int *lda, RealScalar *px, int *incx)
271 {
272 Scalar* a = reinterpret_cast<Scalar*>(pa);
273 Scalar* x = reinterpret_cast<Scalar*>(px);
274 int coeff_rows = *k + 1;
275
276 int info = 0;
277 if(UPLO(*uplo)==INVALID) info = 1;
278 else if(OP(*opa)==INVALID) info = 2;
279 else if(DIAG(*diag)==INVALID) info = 3;
280 else if(*n<0) info = 4;
281 else if(*k<0) info = 5;
282 else if(*lda<coeff_rows) info = 7;
283 else if(*incx==0) info = 9;
284 if(info)
285 return xerbla_(SCALAR_SUFFIX_UP"TBMV ",&info,6);
286
287 if(*n==0)
288 return 0;
289
290 int actual_n = *n;
291
292 Scalar* actual_x = get_compact_vector(x,actual_n,*incx);
293
294 MatrixType mat_coeffs(a,coeff_rows,*n,*lda);
295
296 int ku = UPLO(*uplo)==UPPER ? *k : 0;
297 int kl = UPLO(*uplo)==LOWER ? *k : 0;
298
299 for(int j=0; j<*n; ++j)
300 {
301 int start = std::max(0,j - ku);
302 int end = std::min((*m)-1,j + kl);
303 int len = end - start + 1;
304 int offset = (ku) - j + start;
305
306 if(OP(*trans)==NOTR)
307 vector(actual_y+start,len) += (alpha*actual_x[j]) * mat_coeffs.col(j).segment(offset,len);
308 else if(OP(*trans)==TR)
309 actual_y[j] += alpha * ( mat_coeffs.col(j).segment(offset,len).transpose() * vector(actual_x+start,len) ).value();
310 else
311 actual_y[j] += alpha * ( mat_coeffs.col(j).segment(offset,len).adjoint() * vector(actual_x+start,len) ).value();
312 }
313
314 if(actual_x!=x) delete[] actual_x;
315 if(actual_y!=y) delete[] copy_back(actual_y,y,actual_m,*incy);
316
317 return 0;
318 }
319 #endif
320
321 /** DTBSV solves one of the systems of equations
322 *
323 * A*x = b, or A'*x = b,
324 *
325 * where b and x are n element vectors and A is an n by n unit, or
326 * non-unit, upper or lower triangular band matrix, with ( k + 1 )
327 * diagonals.
328 *
329 * No test for singularity or near-singularity is included in this
330 * routine. Such tests must be performed before calling this routine.
331 */
EIGEN_BLAS_FUNC(tbsv)332 int EIGEN_BLAS_FUNC(tbsv)(char *uplo, char *op, char *diag, int *n, int *k, RealScalar *pa, int *lda, RealScalar *px, int *incx)
333 {
334 typedef void (*functype)(int, int, const Scalar *, int, Scalar *);
335 static functype func[16];
336
337 static bool init = false;
338 if(!init)
339 {
340 for(int k=0; k<16; ++k)
341 func[k] = 0;
342
343 func[NOTR | (UP << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector<int,Upper|0, Scalar,false,Scalar,ColMajor>::run);
344 func[TR | (UP << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector<int,Lower|0, Scalar,false,Scalar,RowMajor>::run);
345 func[ADJ | (UP << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector<int,Lower|0, Scalar,Conj, Scalar,RowMajor>::run);
346
347 func[NOTR | (LO << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector<int,Lower|0, Scalar,false,Scalar,ColMajor>::run);
348 func[TR | (LO << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector<int,Upper|0, Scalar,false,Scalar,RowMajor>::run);
349 func[ADJ | (LO << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector<int,Upper|0, Scalar,Conj, Scalar,RowMajor>::run);
350
351 func[NOTR | (UP << 2) | (UNIT << 3)] = (internal::band_solve_triangular_selector<int,Upper|UnitDiag,Scalar,false,Scalar,ColMajor>::run);
352 func[TR | (UP << 2) | (UNIT << 3)] = (internal::band_solve_triangular_selector<int,Lower|UnitDiag,Scalar,false,Scalar,RowMajor>::run);
353 func[ADJ | (UP << 2) | (UNIT << 3)] = (internal::band_solve_triangular_selector<int,Lower|UnitDiag,Scalar,Conj, Scalar,RowMajor>::run);
354
355 func[NOTR | (LO << 2) | (UNIT << 3)] = (internal::band_solve_triangular_selector<int,Lower|UnitDiag,Scalar,false,Scalar,ColMajor>::run);
356 func[TR | (LO << 2) | (UNIT << 3)] = (internal::band_solve_triangular_selector<int,Upper|UnitDiag,Scalar,false,Scalar,RowMajor>::run);
357 func[ADJ | (LO << 2) | (UNIT << 3)] = (internal::band_solve_triangular_selector<int,Upper|UnitDiag,Scalar,Conj, Scalar,RowMajor>::run);
358
359 init = true;
360 }
361
362 Scalar* a = reinterpret_cast<Scalar*>(pa);
363 Scalar* x = reinterpret_cast<Scalar*>(px);
364 int coeff_rows = *k+1;
365
366 int info = 0;
367 if(UPLO(*uplo)==INVALID) info = 1;
368 else if(OP(*op)==INVALID) info = 2;
369 else if(DIAG(*diag)==INVALID) info = 3;
370 else if(*n<0) info = 4;
371 else if(*k<0) info = 5;
372 else if(*lda<coeff_rows) info = 7;
373 else if(*incx==0) info = 9;
374 if(info)
375 return xerbla_(SCALAR_SUFFIX_UP"TBSV ",&info,6);
376
377 if(*n==0 || (*k==0 && DIAG(*diag)==UNIT))
378 return 0;
379
380 int actual_n = *n;
381
382 Scalar* actual_x = get_compact_vector(x,actual_n,*incx);
383
384 int code = OP(*op) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
385 if(code>=16 || func[code]==0)
386 return 0;
387
388 func[code](*n, *k, a, *lda, actual_x);
389
390 if(actual_x!=x) delete[] copy_back(actual_x,x,actual_n,*incx);
391
392 return 0;
393 }
394
395 /** DTPMV performs one of the matrix-vector operations
396 *
397 * x := A*x, or x := A'*x,
398 *
399 * where x is an n element vector and A is an n by n unit, or non-unit,
400 * upper or lower triangular matrix, supplied in packed form.
401 */
EIGEN_BLAS_FUNC(tpmv)402 int EIGEN_BLAS_FUNC(tpmv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pap, RealScalar *px, int *incx)
403 {
404 typedef void (*functype)(int, const Scalar*, const Scalar*, Scalar*, Scalar);
405 static functype func[16];
406
407 static bool init = false;
408 if(!init)
409 {
410 for(int k=0; k<16; ++k)
411 func[k] = 0;
412
413 func[NOTR | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|0, Scalar,false,Scalar,false,ColMajor>::run);
414 func[TR | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|0, Scalar,false,Scalar,false,RowMajor>::run);
415 func[ADJ | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|0, Scalar,Conj, Scalar,false,RowMajor>::run);
416
417 func[NOTR | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|0, Scalar,false,Scalar,false,ColMajor>::run);
418 func[TR | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|0, Scalar,false,Scalar,false,RowMajor>::run);
419 func[ADJ | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|0, Scalar,Conj, Scalar,false,RowMajor>::run);
420
421 func[NOTR | (UP << 2) | (UNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run);
422 func[TR | (UP << 2) | (UNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run);
423 func[ADJ | (UP << 2) | (UNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run);
424
425 func[NOTR | (LO << 2) | (UNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run);
426 func[TR | (LO << 2) | (UNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run);
427 func[ADJ | (LO << 2) | (UNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run);
428
429 init = true;
430 }
431
432 Scalar* ap = reinterpret_cast<Scalar*>(pap);
433 Scalar* x = reinterpret_cast<Scalar*>(px);
434
435 int info = 0;
436 if(UPLO(*uplo)==INVALID) info = 1;
437 else if(OP(*opa)==INVALID) info = 2;
438 else if(DIAG(*diag)==INVALID) info = 3;
439 else if(*n<0) info = 4;
440 else if(*incx==0) info = 7;
441 if(info)
442 return xerbla_(SCALAR_SUFFIX_UP"TPMV ",&info,6);
443
444 if(*n==0)
445 return 1;
446
447 Scalar* actual_x = get_compact_vector(x,*n,*incx);
448 Matrix<Scalar,Dynamic,1> res(*n);
449 res.setZero();
450
451 int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
452 if(code>=16 || func[code]==0)
453 return 0;
454
455 func[code](*n, ap, actual_x, res.data(), Scalar(1));
456
457 copy_back(res.data(),x,*n,*incx);
458 if(actual_x!=x) delete[] actual_x;
459
460 return 1;
461 }
462
463 /** DTPSV solves one of the systems of equations
464 *
465 * A*x = b, or A'*x = b,
466 *
467 * where b and x are n element vectors and A is an n by n unit, or
468 * non-unit, upper or lower triangular matrix, supplied in packed form.
469 *
470 * No test for singularity or near-singularity is included in this
471 * routine. Such tests must be performed before calling this routine.
472 */
EIGEN_BLAS_FUNC(tpsv)473 int EIGEN_BLAS_FUNC(tpsv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pap, RealScalar *px, int *incx)
474 {
475 typedef void (*functype)(int, const Scalar*, Scalar*);
476 static functype func[16];
477
478 static bool init = false;
479 if(!init)
480 {
481 for(int k=0; k<16; ++k)
482 func[k] = 0;
483
484 func[NOTR | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0, false,ColMajor>::run);
485 func[TR | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0, false,RowMajor>::run);
486 func[ADJ | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0, Conj, RowMajor>::run);
487
488 func[NOTR | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0, false,ColMajor>::run);
489 func[TR | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0, false,RowMajor>::run);
490 func[ADJ | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0, Conj, RowMajor>::run);
491
492 func[NOTR | (UP << 2) | (UNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,ColMajor>::run);
493 func[TR | (UP << 2) | (UNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,RowMajor>::run);
494 func[ADJ | (UP << 2) | (UNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,Conj, RowMajor>::run);
495
496 func[NOTR | (LO << 2) | (UNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,ColMajor>::run);
497 func[TR | (LO << 2) | (UNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,RowMajor>::run);
498 func[ADJ | (LO << 2) | (UNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,Conj, RowMajor>::run);
499
500 init = true;
501 }
502
503 Scalar* ap = reinterpret_cast<Scalar*>(pap);
504 Scalar* x = reinterpret_cast<Scalar*>(px);
505
506 int info = 0;
507 if(UPLO(*uplo)==INVALID) info = 1;
508 else if(OP(*opa)==INVALID) info = 2;
509 else if(DIAG(*diag)==INVALID) info = 3;
510 else if(*n<0) info = 4;
511 else if(*incx==0) info = 7;
512 if(info)
513 return xerbla_(SCALAR_SUFFIX_UP"TPSV ",&info,6);
514
515 Scalar* actual_x = get_compact_vector(x,*n,*incx);
516
517 int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
518 func[code](*n, ap, actual_x);
519
520 if(actual_x!=x) delete[] copy_back(actual_x,x,*n,*incx);
521
522 return 1;
523 }
524
525