1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library
2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra.
3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath//
4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2009-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath//
6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla
7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed
8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath
10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include "common.h"
11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath
12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// y = alpha*A*x + beta*y
13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathint EIGEN_BLAS_FUNC(symv) (char *uplo, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *px, int *incx, RealScalar *pbeta, RealScalar *py, int *incy)
14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{
15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  typedef void (*functype)(int, const Scalar*, int, const Scalar*, int, Scalar*, Scalar);
167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez  static functype func[2];
17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath
18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  static bool init = false;
19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  if(!init)
20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  {
21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath    for(int k=0; k<2; ++k)
22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath      func[k] = 0;
23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath
24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath    func[UP] = (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Upper,false,false>::run);
25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath    func[LO] = (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Lower,false,false>::run);
26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath
27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath    init = true;
28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  }
29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath
30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  Scalar* a = reinterpret_cast<Scalar*>(pa);
31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  Scalar* x = reinterpret_cast<Scalar*>(px);
32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  Scalar* y = reinterpret_cast<Scalar*>(py);
33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  Scalar alpha  = *reinterpret_cast<Scalar*>(palpha);
34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath  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
71int 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
132int 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  */
234int 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  */
285int 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  */
340int 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
371