RavEngine
Loading...
Searching...
No Matches
level2_cplx_impl.h
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
19int 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
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
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
105int 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
149int 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
197int 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
244int 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
295int 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
332int 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}
Eigen::selfadjoint_rank1_update
Definition GeneralMatrixMatrixTriangular.h:16
internal::general_rank1_update
Definition GeneralRank1Update.h:17
internal::packed_rank2_update_selector
Definition Rank2Update.h:38
internal::rank2_update_selector
Definition Rank2Update.h:20
internal::selfadjoint_packed_rank1_update
Definition PackedSelfadjointProduct.h:19