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