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 | // y = alpha*A*x + beta*y
|
---|
13 | int EIGEN_BLAS_FUNC(symv) (char *uplo, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *px, int *incx, RealScalar *pbeta, RealScalar *py, int *incy)
|
---|
14 | {
|
---|
15 | typedef void (*functype)(int, const Scalar*, int, const Scalar*, int, Scalar*, Scalar);
|
---|
16 | static functype func[2];
|
---|
17 |
|
---|
18 | static bool init = false;
|
---|
19 | if(!init)
|
---|
20 | {
|
---|
21 | for(int k=0; k<2; ++k)
|
---|
22 | func[k] = 0;
|
---|
23 |
|
---|
24 | func[UP] = (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Upper,false,false>::run);
|
---|
25 | func[LO] = (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Lower,false,false>::run);
|
---|
26 |
|
---|
27 | init = true;
|
---|
28 | }
|
---|
29 |
|
---|
30 | Scalar* a = reinterpret_cast<Scalar*>(pa);
|
---|
31 | Scalar* x = reinterpret_cast<Scalar*>(px);
|
---|
32 | Scalar* y = reinterpret_cast<Scalar*>(py);
|
---|
33 | Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
|
---|
34 | Scalar beta = *reinterpret_cast<Scalar*>(pbeta);
|
---|
35 |
|
---|
36 | // check arguments
|
---|
37 | int info = 0;
|
---|
38 | if(UPLO(*uplo)==INVALID) info = 1;
|
---|
39 | else if(*n<0) info = 2;
|
---|
40 | else if(*lda<std::max(1,*n)) info = 5;
|
---|
41 | else if(*incx==0) info = 7;
|
---|
42 | else if(*incy==0) info = 10;
|
---|
43 | if(info)
|
---|
44 | return xerbla_(SCALAR_SUFFIX_UP"SYMV ",&info,6);
|
---|
45 |
|
---|
46 | if(*n==0)
|
---|
47 | return 0;
|
---|
48 |
|
---|
49 | Scalar* actual_x = get_compact_vector(x,*n,*incx);
|
---|
50 | Scalar* actual_y = get_compact_vector(y,*n,*incy);
|
---|
51 |
|
---|
52 | if(beta!=Scalar(1))
|
---|
53 | {
|
---|
54 | if(beta==Scalar(0)) vector(actual_y, *n).setZero();
|
---|
55 | else vector(actual_y, *n) *= beta;
|
---|
56 | }
|
---|
57 |
|
---|
58 | int code = UPLO(*uplo);
|
---|
59 | if(code>=2 || func[code]==0)
|
---|
60 | return 0;
|
---|
61 |
|
---|
62 | func[code](*n, a, *lda, actual_x, 1, actual_y, alpha);
|
---|
63 |
|
---|
64 | if(actual_x!=x) delete[] actual_x;
|
---|
65 | if(actual_y!=y) delete[] copy_back(actual_y,y,*n,*incy);
|
---|
66 |
|
---|
67 | return 1;
|
---|
68 | }
|
---|
69 |
|
---|
70 | // C := alpha*x*x' + C
|
---|
71 | int EIGEN_BLAS_FUNC(syr)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *pc, int *ldc)
|
---|
72 | {
|
---|
73 |
|
---|
74 | // typedef void (*functype)(int, const Scalar *, int, Scalar *, int, Scalar);
|
---|
75 | // static functype func[2];
|
---|
76 |
|
---|
77 | // static bool init = false;
|
---|
78 | // if(!init)
|
---|
79 | // {
|
---|
80 | // for(int k=0; k<2; ++k)
|
---|
81 | // func[k] = 0;
|
---|
82 | //
|
---|
83 | // func[UP] = (internal::selfadjoint_product<Scalar,ColMajor,ColMajor,false,UpperTriangular>::run);
|
---|
84 | // func[LO] = (internal::selfadjoint_product<Scalar,ColMajor,ColMajor,false,LowerTriangular>::run);
|
---|
85 |
|
---|
86 | // init = true;
|
---|
87 | // }
|
---|
88 | typedef void (*functype)(int, Scalar*, int, const Scalar*, const Scalar*, const Scalar&);
|
---|
89 | static functype func[2];
|
---|
90 |
|
---|
91 | static bool init = false;
|
---|
92 | if(!init)
|
---|
93 | {
|
---|
94 | for(int k=0; k<2; ++k)
|
---|
95 | func[k] = 0;
|
---|
96 |
|
---|
97 | func[UP] = (selfadjoint_rank1_update<Scalar,int,ColMajor,Upper,false,Conj>::run);
|
---|
98 | func[LO] = (selfadjoint_rank1_update<Scalar,int,ColMajor,Lower,false,Conj>::run);
|
---|
99 |
|
---|
100 | init = true;
|
---|
101 | }
|
---|
102 |
|
---|
103 | Scalar* x = reinterpret_cast<Scalar*>(px);
|
---|
104 | Scalar* c = reinterpret_cast<Scalar*>(pc);
|
---|
105 | Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
|
---|
106 |
|
---|
107 | int info = 0;
|
---|
108 | if(UPLO(*uplo)==INVALID) info = 1;
|
---|
109 | else if(*n<0) info = 2;
|
---|
110 | else if(*incx==0) info = 5;
|
---|
111 | else if(*ldc<std::max(1,*n)) info = 7;
|
---|
112 | if(info)
|
---|
113 | return xerbla_(SCALAR_SUFFIX_UP"SYR ",&info,6);
|
---|
114 |
|
---|
115 | if(*n==0 || alpha==Scalar(0)) return 1;
|
---|
116 |
|
---|
117 | // if the increment is not 1, let's copy it to a temporary vector to enable vectorization
|
---|
118 | Scalar* x_cpy = get_compact_vector(x,*n,*incx);
|
---|
119 |
|
---|
120 | int code = UPLO(*uplo);
|
---|
121 | if(code>=2 || func[code]==0)
|
---|
122 | return 0;
|
---|
123 |
|
---|
124 | func[code](*n, c, *ldc, x_cpy, x_cpy, alpha);
|
---|
125 |
|
---|
126 | if(x_cpy!=x) delete[] x_cpy;
|
---|
127 |
|
---|
128 | return 1;
|
---|
129 | }
|
---|
130 |
|
---|
131 | // C := alpha*x*y' + alpha*y*x' + C
|
---|
132 | int EIGEN_BLAS_FUNC(syr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pc, int *ldc)
|
---|
133 | {
|
---|
134 | // typedef void (*functype)(int, const Scalar *, int, const Scalar *, int, Scalar *, int, Scalar);
|
---|
135 | // static functype func[2];
|
---|
136 | //
|
---|
137 | // static bool init = false;
|
---|
138 | // if(!init)
|
---|
139 | // {
|
---|
140 | // for(int k=0; k<2; ++k)
|
---|
141 | // func[k] = 0;
|
---|
142 | //
|
---|
143 | // func[UP] = (internal::selfadjoint_product<Scalar,ColMajor,ColMajor,false,UpperTriangular>::run);
|
---|
144 | // func[LO] = (internal::selfadjoint_product<Scalar,ColMajor,ColMajor,false,LowerTriangular>::run);
|
---|
145 | //
|
---|
146 | // init = true;
|
---|
147 | // }
|
---|
148 | typedef void (*functype)(int, Scalar*, int, const Scalar*, const Scalar*, Scalar);
|
---|
149 | static functype func[2];
|
---|
150 |
|
---|
151 | static bool init = false;
|
---|
152 | if(!init)
|
---|
153 | {
|
---|
154 | for(int k=0; k<2; ++k)
|
---|
155 | func[k] = 0;
|
---|
156 |
|
---|
157 | func[UP] = (internal::rank2_update_selector<Scalar,int,Upper>::run);
|
---|
158 | func[LO] = (internal::rank2_update_selector<Scalar,int,Lower>::run);
|
---|
159 |
|
---|
160 | init = true;
|
---|
161 | }
|
---|
162 |
|
---|
163 | Scalar* x = reinterpret_cast<Scalar*>(px);
|
---|
164 | Scalar* y = reinterpret_cast<Scalar*>(py);
|
---|
165 | Scalar* c = reinterpret_cast<Scalar*>(pc);
|
---|
166 | Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
|
---|
167 |
|
---|
168 | int info = 0;
|
---|
169 | if(UPLO(*uplo)==INVALID) info = 1;
|
---|
170 | else if(*n<0) info = 2;
|
---|
171 | else if(*incx==0) info = 5;
|
---|
172 | else if(*incy==0) info = 7;
|
---|
173 | else if(*ldc<std::max(1,*n)) info = 9;
|
---|
174 | if(info)
|
---|
175 | return xerbla_(SCALAR_SUFFIX_UP"SYR2 ",&info,6);
|
---|
176 |
|
---|
177 | if(alpha==Scalar(0))
|
---|
178 | return 1;
|
---|
179 |
|
---|
180 | Scalar* x_cpy = get_compact_vector(x,*n,*incx);
|
---|
181 | Scalar* y_cpy = get_compact_vector(y,*n,*incy);
|
---|
182 |
|
---|
183 | int code = UPLO(*uplo);
|
---|
184 | if(code>=2 || func[code]==0)
|
---|
185 | return 0;
|
---|
186 |
|
---|
187 | func[code](*n, c, *ldc, x_cpy, y_cpy, alpha);
|
---|
188 |
|
---|
189 | if(x_cpy!=x) delete[] x_cpy;
|
---|
190 | if(y_cpy!=y) delete[] y_cpy;
|
---|
191 |
|
---|
192 | // int code = UPLO(*uplo);
|
---|
193 | // if(code>=2 || func[code]==0)
|
---|
194 | // return 0;
|
---|
195 |
|
---|
196 | // func[code](*n, a, *inca, b, *incb, c, *ldc, alpha);
|
---|
197 | return 1;
|
---|
198 | }
|
---|
199 |
|
---|
200 | /** DSBMV performs the matrix-vector operation
|
---|
201 | *
|
---|
202 | * y := alpha*A*x + beta*y,
|
---|
203 | *
|
---|
204 | * where alpha and beta are scalars, x and y are n element vectors and
|
---|
205 | * A is an n by n symmetric band matrix, with k super-diagonals.
|
---|
206 | */
|
---|
207 | // int EIGEN_BLAS_FUNC(sbmv)( char *uplo, int *n, int *k, RealScalar *alpha, RealScalar *a, int *lda,
|
---|
208 | // RealScalar *x, int *incx, RealScalar *beta, RealScalar *y, int *incy)
|
---|
209 | // {
|
---|
210 | // return 1;
|
---|
211 | // }
|
---|
212 |
|
---|
213 |
|
---|
214 | /** DSPMV performs the matrix-vector operation
|
---|
215 | *
|
---|
216 | * y := alpha*A*x + beta*y,
|
---|
217 | *
|
---|
218 | * where alpha and beta are scalars, x and y are n element vectors and
|
---|
219 | * A is an n by n symmetric matrix, supplied in packed form.
|
---|
220 | *
|
---|
221 | */
|
---|
222 | // int EIGEN_BLAS_FUNC(spmv)(char *uplo, int *n, RealScalar *alpha, RealScalar *ap, RealScalar *x, int *incx, RealScalar *beta, RealScalar *y, int *incy)
|
---|
223 | // {
|
---|
224 | // return 1;
|
---|
225 | // }
|
---|
226 |
|
---|
227 | /** DSPR performs the symmetric rank 1 operation
|
---|
228 | *
|
---|
229 | * A := alpha*x*x' + A,
|
---|
230 | *
|
---|
231 | * where alpha is a real scalar, x is an n element vector and A is an
|
---|
232 | * n by n symmetric matrix, supplied in packed form.
|
---|
233 | */
|
---|
234 | int EIGEN_BLAS_FUNC(spr)(char *uplo, int *n, Scalar *palpha, Scalar *px, int *incx, Scalar *pap)
|
---|
235 | {
|
---|
236 | typedef void (*functype)(int, Scalar*, const Scalar*, Scalar);
|
---|
237 | static functype func[2];
|
---|
238 |
|
---|
239 | static bool init = false;
|
---|
240 | if(!init)
|
---|
241 | {
|
---|
242 | for(int k=0; k<2; ++k)
|
---|
243 | func[k] = 0;
|
---|
244 |
|
---|
245 | func[UP] = (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Upper,false,false>::run);
|
---|
246 | func[LO] = (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Lower,false,false>::run);
|
---|
247 |
|
---|
248 | init = true;
|
---|
249 | }
|
---|
250 |
|
---|
251 | Scalar* x = reinterpret_cast<Scalar*>(px);
|
---|
252 | Scalar* ap = reinterpret_cast<Scalar*>(pap);
|
---|
253 | Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
|
---|
254 |
|
---|
255 | int info = 0;
|
---|
256 | if(UPLO(*uplo)==INVALID) info = 1;
|
---|
257 | else if(*n<0) info = 2;
|
---|
258 | else if(*incx==0) info = 5;
|
---|
259 | if(info)
|
---|
260 | return xerbla_(SCALAR_SUFFIX_UP"SPR ",&info,6);
|
---|
261 |
|
---|
262 | if(alpha==Scalar(0))
|
---|
263 | return 1;
|
---|
264 |
|
---|
265 | Scalar* x_cpy = get_compact_vector(x, *n, *incx);
|
---|
266 |
|
---|
267 | int code = UPLO(*uplo);
|
---|
268 | if(code>=2 || func[code]==0)
|
---|
269 | return 0;
|
---|
270 |
|
---|
271 | func[code](*n, ap, x_cpy, alpha);
|
---|
272 |
|
---|
273 | if(x_cpy!=x) delete[] x_cpy;
|
---|
274 |
|
---|
275 | return 1;
|
---|
276 | }
|
---|
277 |
|
---|
278 | /** DSPR2 performs the symmetric rank 2 operation
|
---|
279 | *
|
---|
280 | * A := alpha*x*y' + alpha*y*x' + A,
|
---|
281 | *
|
---|
282 | * where alpha is a scalar, x and y are n element vectors and A is an
|
---|
283 | * n by n symmetric matrix, supplied in packed form.
|
---|
284 | */
|
---|
285 | int EIGEN_BLAS_FUNC(spr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pap)
|
---|
286 | {
|
---|
287 | typedef void (*functype)(int, Scalar*, const Scalar*, const Scalar*, Scalar);
|
---|
288 | static functype func[2];
|
---|
289 |
|
---|
290 | static bool init = false;
|
---|
291 | if(!init)
|
---|
292 | {
|
---|
293 | for(int k=0; k<2; ++k)
|
---|
294 | func[k] = 0;
|
---|
295 |
|
---|
296 | func[UP] = (internal::packed_rank2_update_selector<Scalar,int,Upper>::run);
|
---|
297 | func[LO] = (internal::packed_rank2_update_selector<Scalar,int,Lower>::run);
|
---|
298 |
|
---|
299 | init = true;
|
---|
300 | }
|
---|
301 |
|
---|
302 | Scalar* x = reinterpret_cast<Scalar*>(px);
|
---|
303 | Scalar* y = reinterpret_cast<Scalar*>(py);
|
---|
304 | Scalar* ap = reinterpret_cast<Scalar*>(pap);
|
---|
305 | Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
|
---|
306 |
|
---|
307 | int info = 0;
|
---|
308 | if(UPLO(*uplo)==INVALID) info = 1;
|
---|
309 | else if(*n<0) info = 2;
|
---|
310 | else if(*incx==0) info = 5;
|
---|
311 | else if(*incy==0) info = 7;
|
---|
312 | if(info)
|
---|
313 | return xerbla_(SCALAR_SUFFIX_UP"SPR2 ",&info,6);
|
---|
314 |
|
---|
315 | if(alpha==Scalar(0))
|
---|
316 | return 1;
|
---|
317 |
|
---|
318 | Scalar* x_cpy = get_compact_vector(x, *n, *incx);
|
---|
319 | Scalar* y_cpy = get_compact_vector(y, *n, *incy);
|
---|
320 |
|
---|
321 | int code = UPLO(*uplo);
|
---|
322 | if(code>=2 || func[code]==0)
|
---|
323 | return 0;
|
---|
324 |
|
---|
325 | func[code](*n, ap, x_cpy, y_cpy, alpha);
|
---|
326 |
|
---|
327 | if(x_cpy!=x) delete[] x_cpy;
|
---|
328 | if(y_cpy!=y) delete[] y_cpy;
|
---|
329 |
|
---|
330 | return 1;
|
---|
331 | }
|
---|
332 |
|
---|
333 | /** DGER performs the rank 1 operation
|
---|
334 | *
|
---|
335 | * A := alpha*x*y' + A,
|
---|
336 | *
|
---|
337 | * where alpha is a scalar, x is an m element vector, y is an n element
|
---|
338 | * vector and A is an m by n matrix.
|
---|
339 | */
|
---|
340 | int EIGEN_BLAS_FUNC(ger)(int *m, int *n, Scalar *palpha, Scalar *px, int *incx, Scalar *py, int *incy, Scalar *pa, int *lda)
|
---|
341 | {
|
---|
342 | Scalar* x = reinterpret_cast<Scalar*>(px);
|
---|
343 | Scalar* y = reinterpret_cast<Scalar*>(py);
|
---|
344 | Scalar* a = reinterpret_cast<Scalar*>(pa);
|
---|
345 | Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
|
---|
346 |
|
---|
347 | int info = 0;
|
---|
348 | if(*m<0) info = 1;
|
---|
349 | else if(*n<0) info = 2;
|
---|
350 | else if(*incx==0) info = 5;
|
---|
351 | else if(*incy==0) info = 7;
|
---|
352 | else if(*lda<std::max(1,*m)) info = 9;
|
---|
353 | if(info)
|
---|
354 | return xerbla_(SCALAR_SUFFIX_UP"GER ",&info,6);
|
---|
355 |
|
---|
356 | if(alpha==Scalar(0))
|
---|
357 | return 1;
|
---|
358 |
|
---|
359 | Scalar* x_cpy = get_compact_vector(x,*m,*incx);
|
---|
360 | Scalar* y_cpy = get_compact_vector(y,*n,*incy);
|
---|
361 |
|
---|
362 | internal::general_rank1_update<Scalar,int,ColMajor,false,false>::run(*m, *n, a, *lda, x_cpy, y_cpy, alpha);
|
---|
363 |
|
---|
364 | if(x_cpy!=x) delete[] x_cpy;
|
---|
365 | if(y_cpy!=y) delete[] y_cpy;
|
---|
366 |
|
---|
367 | return 1;
|
---|
368 | }
|
---|
369 |
|
---|
370 |
|
---|