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) (char *uplo, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *px, int *incx, RealScalar *pbeta, RealScalar *py, int *incy) 14 { 15 typedef void (*functype)(int, const Scalar*, int, const Scalar*, int, Scalar*, Scalar); 16 static functype func[2]; 17 18 static bool init = false; 19 if(!init) 20 { 21 for(int k=0; k<2; ++k) 22 func[k] = 0; 23 24 func[UP] = (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Upper,false,false>::run); 25 func[LO] = (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Lower,false,false>::run); 26 27 init = true; 28 } 29 30 Scalar* a = reinterpret_cast<Scalar*>(pa); 31 Scalar* x = reinterpret_cast<Scalar*>(px); 32 Scalar* y = reinterpret_cast<Scalar*>(py); 33 Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 34 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 71 int 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 132 int 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 */ 234 int 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 */ 285 int 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 */ 340 int 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
