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