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 (at) 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 13 int EIGEN_BLAS_FUNC(symv) (const char *uplo, const int *n, const RealScalar *palpha, const RealScalar *pa, const int *lda, 14 const RealScalar *px, const int *incx, const RealScalar *pbeta, RealScalar *py, const int *incy) 15 { 16 typedef void (*functype)(int, const Scalar*, int, const Scalar*, Scalar*, Scalar); 17 static const functype func[2] = { 18 // array index: UP 19 (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Upper,false,false>::run), 20 // array index: LO 21 (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Lower,false,false>::run), 22 }; 23 24 const Scalar* a = reinterpret_cast<const Scalar*>(pa); 25 const Scalar* x = reinterpret_cast<const Scalar*>(px); 26 Scalar* y = reinterpret_cast<Scalar*>(py); 27 Scalar alpha = *reinterpret_cast<const Scalar*>(palpha); 28 Scalar beta = *reinterpret_cast<const Scalar*>(pbeta); 29 30 // check arguments 31 int info = 0; 32 if(UPLO(*uplo)==INVALID) info = 1; 33 else if(*n<0) info = 2; 34 else if(*lda<std::max(1,*n)) info = 5; 35 else if(*incx==0) info = 7; 36 else if(*incy==0) info = 10; 37 if(info) 38 return xerbla_(SCALAR_SUFFIX_UP"SYMV ",&info,6); 39 40 if(*n==0) 41 return 0; 42 43 const Scalar* actual_x = get_compact_vector(x,*n,*incx); 44 Scalar* actual_y = get_compact_vector(y,*n,*incy); 45 46 if(beta!=Scalar(1)) 47 { 48 if(beta==Scalar(0)) make_vector(actual_y, *n).setZero(); 49 else make_vector(actual_y, *n) *= beta; 50 } 51 52 int code = UPLO(*uplo); 53 if(code>=2 || func[code]==0) 54 return 0; 55 56 func[code](*n, a, *lda, actual_x, actual_y, alpha); 57 58 if(actual_x!=x) delete[] actual_x; 59 if(actual_y!=y) delete[] copy_back(actual_y,y,*n,*incy); 60 61 return 1; 62 } 63 64 // C := alpha*x*x' + C 65 int EIGEN_BLAS_FUNC(syr)(const char *uplo, const int *n, const RealScalar *palpha, const RealScalar *px, const int *incx, RealScalar *pc, const int *ldc) 66 { 67 68 typedef void (*functype)(int, Scalar*, int, const Scalar*, const Scalar*, const Scalar&); 69 static const functype func[2] = { 70 // array index: UP 71 (selfadjoint_rank1_update<Scalar,int,ColMajor,Upper,false,Conj>::run), 72 // array index: LO 73 (selfadjoint_rank1_update<Scalar,int,ColMajor,Lower,false,Conj>::run), 74 }; 75 76 const Scalar* x = reinterpret_cast<const Scalar*>(px); 77 Scalar* c = reinterpret_cast<Scalar*>(pc); 78 Scalar alpha = *reinterpret_cast<const Scalar*>(palpha); 79 80 int info = 0; 81 if(UPLO(*uplo)==INVALID) info = 1; 82 else if(*n<0) info = 2; 83 else if(*incx==0) info = 5; 84 else if(*ldc<std::max(1,*n)) info = 7; 85 if(info) 86 return xerbla_(SCALAR_SUFFIX_UP"SYR ",&info,6); 87 88 if(*n==0 || alpha==Scalar(0)) return 1; 89 90 // if the increment is not 1, let's copy it to a temporary vector to enable vectorization 91 const Scalar* x_cpy = get_compact_vector(x,*n,*incx); 92 93 int code = UPLO(*uplo); 94 if(code>=2 || func[code]==0) 95 return 0; 96 97 func[code](*n, c, *ldc, x_cpy, x_cpy, alpha); 98 99 if(x_cpy!=x) delete[] x_cpy; 100 101 return 1; 102 } 103 104 // C := alpha*x*y' + alpha*y*x' + C 105 int EIGEN_BLAS_FUNC(syr2)(const char *uplo, const int *n, const RealScalar *palpha, const RealScalar *px, const int *incx, const RealScalar *py, const int *incy, RealScalar *pc, const int *ldc) 106 { 107 typedef void (*functype)(int, Scalar*, int, const Scalar*, const Scalar*, Scalar); 108 static const functype func[2] = { 109 // array index: UP 110 (internal::rank2_update_selector<Scalar,int,Upper>::run), 111 // array index: LO 112 (internal::rank2_update_selector<Scalar,int,Lower>::run), 113 }; 114 115 const Scalar* x = reinterpret_cast<const Scalar*>(px); 116 const Scalar* y = reinterpret_cast<const Scalar*>(py); 117 Scalar* c = reinterpret_cast<Scalar*>(pc); 118 Scalar alpha = *reinterpret_cast<const Scalar*>(palpha); 119 120 int info = 0; 121 if(UPLO(*uplo)==INVALID) info = 1; 122 else if(*n<0) info = 2; 123 else if(*incx==0) info = 5; 124 else if(*incy==0) info = 7; 125 else if(*ldc<std::max(1,*n)) info = 9; 126 if(info) 127 return xerbla_(SCALAR_SUFFIX_UP"SYR2 ",&info,6); 128 129 if(alpha==Scalar(0)) 130 return 1; 131 132 const Scalar* x_cpy = get_compact_vector(x,*n,*incx); 133 const Scalar* y_cpy = get_compact_vector(y,*n,*incy); 134 135 int code = UPLO(*uplo); 136 if(code>=2 || func[code]==0) 137 return 0; 138 139 func[code](*n, c, *ldc, x_cpy, y_cpy, alpha); 140 141 if(x_cpy!=x) delete[] x_cpy; 142 if(y_cpy!=y) delete[] y_cpy; 143 144 // int code = UPLO(*uplo); 145 // if(code>=2 || func[code]==0) 146 // return 0; 147 148 // func[code](*n, a, *inca, b, *incb, c, *ldc, alpha); 149 return 1; 150 } 151 152 /** DSBMV performs the matrix-vector operation 153 * 154 * y := alpha*A*x + beta*y, 155 * 156 * where alpha and beta are scalars, x and y are n element vectors and 157 * A is an n by n symmetric band matrix, with k super-diagonals. 158 */ 159 // int EIGEN_BLAS_FUNC(sbmv)( char *uplo, int *n, int *k, RealScalar *alpha, RealScalar *a, int *lda, 160 // RealScalar *x, int *incx, RealScalar *beta, RealScalar *y, int *incy) 161 // { 162 // return 1; 163 // } 164 165 166 /** DSPMV performs the matrix-vector operation 167 * 168 * y := alpha*A*x + beta*y, 169 * 170 * where alpha and beta are scalars, x and y are n element vectors and 171 * A is an n by n symmetric matrix, supplied in packed form. 172 * 173 */ 174 // int EIGEN_BLAS_FUNC(spmv)(char *uplo, int *n, RealScalar *alpha, RealScalar *ap, RealScalar *x, int *incx, RealScalar *beta, RealScalar *y, int *incy) 175 // { 176 // return 1; 177 // } 178 179 /** DSPR performs the symmetric rank 1 operation 180 * 181 * A := alpha*x*x' + A, 182 * 183 * where alpha is a real scalar, x is an n element vector and A is an 184 * n by n symmetric matrix, supplied in packed form. 185 */ 186 int EIGEN_BLAS_FUNC(spr)(char *uplo, int *n, Scalar *palpha, Scalar *px, int *incx, Scalar *pap) 187 { 188 typedef void (*functype)(int, Scalar*, const Scalar*, Scalar); 189 static const functype func[2] = { 190 // array index: UP 191 (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Upper,false,false>::run), 192 // array index: LO 193 (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Lower,false,false>::run), 194 }; 195 196 Scalar* x = reinterpret_cast<Scalar*>(px); 197 Scalar* ap = reinterpret_cast<Scalar*>(pap); 198 Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 199 200 int info = 0; 201 if(UPLO(*uplo)==INVALID) info = 1; 202 else if(*n<0) info = 2; 203 else if(*incx==0) info = 5; 204 if(info) 205 return xerbla_(SCALAR_SUFFIX_UP"SPR ",&info,6); 206 207 if(alpha==Scalar(0)) 208 return 1; 209 210 Scalar* x_cpy = get_compact_vector(x, *n, *incx); 211 212 int code = UPLO(*uplo); 213 if(code>=2 || func[code]==0) 214 return 0; 215 216 func[code](*n, ap, x_cpy, alpha); 217 218 if(x_cpy!=x) delete[] x_cpy; 219 220 return 1; 221 } 222 223 /** DSPR2 performs the symmetric rank 2 operation 224 * 225 * A := alpha*x*y' + alpha*y*x' + A, 226 * 227 * where alpha is a scalar, x and y are n element vectors and A is an 228 * n by n symmetric matrix, supplied in packed form. 229 */ 230 int EIGEN_BLAS_FUNC(spr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pap) 231 { 232 typedef void (*functype)(int, Scalar*, const Scalar*, const Scalar*, Scalar); 233 static const functype func[2] = { 234 // array index: UP 235 (internal::packed_rank2_update_selector<Scalar,int,Upper>::run), 236 // array index: LO 237 (internal::packed_rank2_update_selector<Scalar,int,Lower>::run), 238 }; 239 240 Scalar* x = reinterpret_cast<Scalar*>(px); 241 Scalar* y = reinterpret_cast<Scalar*>(py); 242 Scalar* ap = reinterpret_cast<Scalar*>(pap); 243 Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 244 245 int info = 0; 246 if(UPLO(*uplo)==INVALID) info = 1; 247 else if(*n<0) info = 2; 248 else if(*incx==0) info = 5; 249 else if(*incy==0) info = 7; 250 if(info) 251 return xerbla_(SCALAR_SUFFIX_UP"SPR2 ",&info,6); 252 253 if(alpha==Scalar(0)) 254 return 1; 255 256 Scalar* x_cpy = get_compact_vector(x, *n, *incx); 257 Scalar* y_cpy = get_compact_vector(y, *n, *incy); 258 259 int code = UPLO(*uplo); 260 if(code>=2 || func[code]==0) 261 return 0; 262 263 func[code](*n, ap, x_cpy, y_cpy, alpha); 264 265 if(x_cpy!=x) delete[] x_cpy; 266 if(y_cpy!=y) delete[] y_cpy; 267 268 return 1; 269 } 270 271 /** DGER performs the rank 1 operation 272 * 273 * A := alpha*x*y' + A, 274 * 275 * where alpha is a scalar, x is an m element vector, y is an n element 276 * vector and A is an m by n matrix. 277 */ 278 int EIGEN_BLAS_FUNC(ger)(int *m, int *n, Scalar *palpha, Scalar *px, int *incx, Scalar *py, int *incy, Scalar *pa, int *lda) 279 { 280 Scalar* x = reinterpret_cast<Scalar*>(px); 281 Scalar* y = reinterpret_cast<Scalar*>(py); 282 Scalar* a = reinterpret_cast<Scalar*>(pa); 283 Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 284 285 int info = 0; 286 if(*m<0) info = 1; 287 else if(*n<0) info = 2; 288 else if(*incx==0) info = 5; 289 else if(*incy==0) info = 7; 290 else if(*lda<std::max(1,*m)) info = 9; 291 if(info) 292 return xerbla_(SCALAR_SUFFIX_UP"GER ",&info,6); 293 294 if(alpha==Scalar(0)) 295 return 1; 296 297 Scalar* x_cpy = get_compact_vector(x,*m,*incx); 298 Scalar* y_cpy = get_compact_vector(y,*n,*incy); 299 300 internal::general_rank1_update<Scalar,int,ColMajor,false,false>::run(*m, *n, a, *lda, x_cpy, y_cpy, alpha); 301 302 if(x_cpy!=x) delete[] x_cpy; 303 if(y_cpy!=y) delete[] y_cpy; 304 305 return 1; 306 } 307
