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