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 */
EIGEN_BLAS_FUNC(hemv)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 */
EIGEN_BLAS_FUNC(hpr)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 */
EIGEN_BLAS_FUNC(hpr2)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 */
EIGEN_BLAS_FUNC(her)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 */
EIGEN_BLAS_FUNC(her2)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 */
EIGEN_BLAS_FUNC(geru)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 */
EIGEN_BLAS_FUNC(gerc)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 }
395