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// y = alpha*A*x + beta*y
13int EIGEN_BLAS_FUNC(symv) (char *uplo, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *px, int *incx, RealScalar *pbeta, RealScalar *py, int *incy)
14{
15  Scalar* a = reinterpret_cast<Scalar*>(pa);
16  Scalar* x = reinterpret_cast<Scalar*>(px);
17  Scalar* y = reinterpret_cast<Scalar*>(py);
18  Scalar alpha  = *reinterpret_cast<Scalar*>(palpha);
19  Scalar beta   = *reinterpret_cast<Scalar*>(pbeta);
20
21  // check arguments
22  int info = 0;
23  if(UPLO(*uplo)==INVALID)        info = 1;
24  else if(*n<0)                   info = 2;
25  else if(*lda<std::max(1,*n))    info = 5;
26  else if(*incx==0)               info = 7;
27  else if(*incy==0)               info = 10;
28  if(info)
29    return xerbla_(SCALAR_SUFFIX_UP"SYMV ",&info,6);
30
31  if(*n==0)
32    return 0;
33
34  Scalar* actual_x = get_compact_vector(x,*n,*incx);
35  Scalar* actual_y = get_compact_vector(y,*n,*incy);
36
37  if(beta!=Scalar(1))
38  {
39    if(beta==Scalar(0)) vector(actual_y, *n).setZero();
40    else                vector(actual_y, *n) *= beta;
41  }
42
43  // TODO performs a direct call to the underlying implementation function
44       if(UPLO(*uplo)==UP) vector(actual_y,*n).noalias() += matrix(a,*n,*n,*lda).selfadjointView<Upper>() * (alpha * vector(actual_x,*n));
45  else if(UPLO(*uplo)==LO) vector(actual_y,*n).noalias() += matrix(a,*n,*n,*lda).selfadjointView<Lower>() * (alpha * vector(actual_x,*n));
46
47  if(actual_x!=x) delete[] actual_x;
48  if(actual_y!=y) delete[] copy_back(actual_y,y,*n,*incy);
49
50  return 1;
51}
52
53// C := alpha*x*x' + C
54int EIGEN_BLAS_FUNC(syr)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *pc, int *ldc)
55{
56
57//   typedef void (*functype)(int, const Scalar *, int, Scalar *, int, Scalar);
58//   static functype func[2];
59
60//   static bool init = false;
61//   if(!init)
62//   {
63//     for(int k=0; k<2; ++k)
64//       func[k] = 0;
65//
66//     func[UP] = (internal::selfadjoint_product<Scalar,ColMajor,ColMajor,false,UpperTriangular>::run);
67//     func[LO] = (internal::selfadjoint_product<Scalar,ColMajor,ColMajor,false,LowerTriangular>::run);
68
69//     init = true;
70//   }
71
72  Scalar* x = reinterpret_cast<Scalar*>(px);
73  Scalar* c = reinterpret_cast<Scalar*>(pc);
74  Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
75
76  int info = 0;
77  if(UPLO(*uplo)==INVALID)                                            info = 1;
78  else if(*n<0)                                                       info = 2;
79  else if(*incx==0)                                                   info = 5;
80  else if(*ldc<std::max(1,*n))                                        info = 7;
81  if(info)
82    return xerbla_(SCALAR_SUFFIX_UP"SYR  ",&info,6);
83
84  if(*n==0 || alpha==Scalar(0)) return 1;
85
86  // if the increment is not 1, let's copy it to a temporary vector to enable vectorization
87  Scalar* x_cpy = get_compact_vector(x,*n,*incx);
88
89  Matrix<Scalar,Dynamic,Dynamic> m2(matrix(c,*n,*n,*ldc));
90
91  // TODO check why this is not accurate enough for lapack tests
92//   if(UPLO(*uplo)==LO)       matrix(c,*n,*n,*ldc).selfadjointView<Lower>().rankUpdate(vector(x_cpy,*n), alpha);
93//   else if(UPLO(*uplo)==UP)  matrix(c,*n,*n,*ldc).selfadjointView<Upper>().rankUpdate(vector(x_cpy,*n), alpha);
94
95  if(UPLO(*uplo)==LO)
96    for(int j=0;j<*n;++j)
97      matrix(c,*n,*n,*ldc).col(j).tail(*n-j) += alpha * x_cpy[j] * vector(x_cpy+j,*n-j);
98  else
99    for(int j=0;j<*n;++j)
100      matrix(c,*n,*n,*ldc).col(j).head(j+1) += alpha * x_cpy[j] * vector(x_cpy,j+1);
101
102  if(x_cpy!=x)  delete[] x_cpy;
103
104  return 1;
105}
106
107// C := alpha*x*y' + alpha*y*x' + C
108int EIGEN_BLAS_FUNC(syr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pc, int *ldc)
109{
110//   typedef void (*functype)(int, const Scalar *, int, const Scalar *, int, Scalar *, int, Scalar);
111//   static functype func[2];
112//
113//   static bool init = false;
114//   if(!init)
115//   {
116//     for(int k=0; k<2; ++k)
117//       func[k] = 0;
118//
119//     func[UP] = (internal::selfadjoint_product<Scalar,ColMajor,ColMajor,false,UpperTriangular>::run);
120//     func[LO] = (internal::selfadjoint_product<Scalar,ColMajor,ColMajor,false,LowerTriangular>::run);
121//
122//     init = true;
123//   }
124
125  Scalar* x = reinterpret_cast<Scalar*>(px);
126  Scalar* y = reinterpret_cast<Scalar*>(py);
127  Scalar* c = reinterpret_cast<Scalar*>(pc);
128  Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
129
130  int info = 0;
131  if(UPLO(*uplo)==INVALID)                                            info = 1;
132  else if(*n<0)                                                       info = 2;
133  else if(*incx==0)                                                   info = 5;
134  else if(*incy==0)                                                   info = 7;
135  else if(*ldc<std::max(1,*n))                                        info = 9;
136  if(info)
137    return xerbla_(SCALAR_SUFFIX_UP"SYR2 ",&info,6);
138
139  if(alpha==Scalar(0))
140    return 1;
141
142  Scalar* x_cpy = get_compact_vector(x,*n,*incx);
143  Scalar* y_cpy = get_compact_vector(y,*n,*incy);
144
145  // TODO perform direct calls to underlying implementation
146  if(UPLO(*uplo)==LO)       matrix(c,*n,*n,*ldc).selfadjointView<Lower>().rankUpdate(vector(x_cpy,*n), vector(y_cpy,*n), alpha);
147  else if(UPLO(*uplo)==UP)  matrix(c,*n,*n,*ldc).selfadjointView<Upper>().rankUpdate(vector(x_cpy,*n), vector(y_cpy,*n), alpha);
148
149  if(x_cpy!=x)  delete[] x_cpy;
150  if(y_cpy!=y)  delete[] y_cpy;
151
152//   int code = UPLO(*uplo);
153//   if(code>=2 || func[code]==0)
154//     return 0;
155
156//   func[code](*n, a, *inca, b, *incb, c, *ldc, alpha);
157  return 1;
158}
159
160/**  DSBMV  performs the matrix-vector  operation
161  *
162  *     y := alpha*A*x + beta*y,
163  *
164  *  where alpha and beta are scalars, x and y are n element vectors and
165  *  A is an n by n symmetric band matrix, with k super-diagonals.
166  */
167// int EIGEN_BLAS_FUNC(sbmv)( char *uplo, int *n, int *k, RealScalar *alpha, RealScalar *a, int *lda,
168//                            RealScalar *x, int *incx, RealScalar *beta, RealScalar *y, int *incy)
169// {
170//   return 1;
171// }
172
173
174/**  DSPMV  performs the matrix-vector operation
175  *
176  *     y := alpha*A*x + beta*y,
177  *
178  *  where alpha and beta are scalars, x and y are n element vectors and
179  *  A is an n by n symmetric matrix, supplied in packed form.
180  *
181  */
182// int EIGEN_BLAS_FUNC(spmv)(char *uplo, int *n, RealScalar *alpha, RealScalar *ap, RealScalar *x, int *incx, RealScalar *beta, RealScalar *y, int *incy)
183// {
184//   return 1;
185// }
186
187/**  DSPR    performs the symmetric rank 1 operation
188  *
189  *     A := alpha*x*x' + A,
190  *
191  *  where alpha is a real scalar, x is an n element vector and A is an
192  *  n by n symmetric matrix, supplied in packed form.
193  */
194// int EIGEN_BLAS_FUNC(spr)(char *uplo, int *n, Scalar *alpha, Scalar *x, int *incx, Scalar *ap)
195// {
196//   return 1;
197// }
198
199/**  DSPR2  performs the symmetric rank 2 operation
200  *
201  *     A := alpha*x*y' + alpha*y*x' + A,
202  *
203  *  where alpha is a scalar, x and y are n element vectors and A is an
204  *  n by n symmetric matrix, supplied in packed form.
205  */
206// int EIGEN_BLAS_FUNC(spr2)(char *uplo, int *n, RealScalar *alpha, RealScalar *x, int *incx, RealScalar *y, int *incy, RealScalar *ap)
207// {
208//   return 1;
209// }
210
211