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