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