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 /** ZHEMV performs the matrix-vector operation 13 * 14 * y := alpha*A*x + beta*y, 15 * 16 * where alpha and beta are scalars, x and y are n element vectors and 17 * A is an n by n hermitian matrix. 18 */ 19 int EIGEN_BLAS_FUNC(hemv)(char *uplo, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *px, int *incx, RealScalar *pbeta, RealScalar *py, int *incy) 20 { 21 typedef void (*functype)(int, const Scalar*, int, const Scalar*, int, Scalar*, Scalar); 22 static functype func[2]; 23 24 static bool init = false; 25 if(!init) 26 { 27 for(int k=0; k<2; ++k) 28 func[k] = 0; 29 30 func[UP] = (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Upper,false,false>::run); 31 func[LO] = (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Lower,false,false>::run); 32 33 init = true; 34 } 35 36 Scalar* a = reinterpret_cast<Scalar*>(pa); 37 Scalar* x = reinterpret_cast<Scalar*>(px); 38 Scalar* y = reinterpret_cast<Scalar*>(py); 39 Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 40 Scalar beta = *reinterpret_cast<Scalar*>(pbeta); 41 42 // check arguments 43 int info = 0; 44 if(UPLO(*uplo)==INVALID) info = 1; 45 else if(*n<0) info = 2; 46 else if(*lda<std::max(1,*n)) info = 5; 47 else if(*incx==0) info = 7; 48 else if(*incy==0) info = 10; 49 if(info) 50 return xerbla_(SCALAR_SUFFIX_UP"HEMV ",&info,6); 51 52 if(*n==0) 53 return 1; 54 55 Scalar* actual_x = get_compact_vector(x,*n,*incx); 56 Scalar* actual_y = get_compact_vector(y,*n,*incy); 57 58 if(beta!=Scalar(1)) 59 { 60 if(beta==Scalar(0)) vector(actual_y, *n).setZero(); 61 else vector(actual_y, *n) *= beta; 62 } 63 64 if(alpha!=Scalar(0)) 65 { 66 int code = UPLO(*uplo); 67 if(code>=2 || func[code]==0) 68 return 0; 69 70 func[code](*n, a, *lda, actual_x, 1, actual_y, alpha); 71 } 72 73 if(actual_x!=x) delete[] actual_x; 74 if(actual_y!=y) delete[] copy_back(actual_y,y,*n,*incy); 75 76 return 1; 77 } 78 79 /** ZHBMV performs the matrix-vector operation 80 * 81 * y := alpha*A*x + beta*y, 82 * 83 * where alpha and beta are scalars, x and y are n element vectors and 84 * A is an n by n hermitian band matrix, with k super-diagonals. 85 */ 86 // int EIGEN_BLAS_FUNC(hbmv)(char *uplo, int *n, int *k, RealScalar *alpha, RealScalar *a, int *lda, 87 // RealScalar *x, int *incx, RealScalar *beta, RealScalar *y, int *incy) 88 // { 89 // return 1; 90 // } 91 92 /** ZHPMV performs the matrix-vector operation 93 * 94 * y := alpha*A*x + beta*y, 95 * 96 * where alpha and beta are scalars, x and y are n element vectors and 97 * A is an n by n hermitian matrix, supplied in packed form. 98 */ 99 // int EIGEN_BLAS_FUNC(hpmv)(char *uplo, int *n, RealScalar *alpha, RealScalar *ap, RealScalar *x, int *incx, RealScalar *beta, RealScalar *y, int *incy) 100 // { 101 // return 1; 102 // } 103 104 /** ZHPR performs the hermitian rank 1 operation 105 * 106 * A := alpha*x*conjg( x' ) + A, 107 * 108 * where alpha is a real scalar, x is an n element vector and A is an 109 * n by n hermitian matrix, supplied in packed form. 110 */ 111 int EIGEN_BLAS_FUNC(hpr)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *pap) 112 { 113 typedef void (*functype)(int, Scalar*, const Scalar*, RealScalar); 114 static functype func[2]; 115 116 static bool init = false; 117 if(!init) 118 { 119 for(int k=0; k<2; ++k) 120 func[k] = 0; 121 122 func[UP] = (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Upper,false,Conj>::run); 123 func[LO] = (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Lower,false,Conj>::run); 124 125 init = true; 126 } 127 128 Scalar* x = reinterpret_cast<Scalar*>(px); 129 Scalar* ap = reinterpret_cast<Scalar*>(pap); 130 RealScalar alpha = *palpha; 131 132 int info = 0; 133 if(UPLO(*uplo)==INVALID) info = 1; 134 else if(*n<0) info = 2; 135 else if(*incx==0) info = 5; 136 if(info) 137 return xerbla_(SCALAR_SUFFIX_UP"HPR ",&info,6); 138 139 if(alpha==Scalar(0)) 140 return 1; 141 142 Scalar* x_cpy = get_compact_vector(x, *n, *incx); 143 144 int code = UPLO(*uplo); 145 if(code>=2 || func[code]==0) 146 return 0; 147 148 func[code](*n, ap, x_cpy, alpha); 149 150 if(x_cpy!=x) delete[] x_cpy; 151 152 return 1; 153 } 154 155 /** ZHPR2 performs the hermitian rank 2 operation 156 * 157 * A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A, 158 * 159 * where alpha is a scalar, x and y are n element vectors and A is an 160 * n by n hermitian matrix, supplied in packed form. 161 */ 162 int EIGEN_BLAS_FUNC(hpr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pap) 163 { 164 typedef void (*functype)(int, Scalar*, const Scalar*, const Scalar*, Scalar); 165 static functype func[2]; 166 167 static bool init = false; 168 if(!init) 169 { 170 for(int k=0; k<2; ++k) 171 func[k] = 0; 172 173 func[UP] = (internal::packed_rank2_update_selector<Scalar,int,Upper>::run); 174 func[LO] = (internal::packed_rank2_update_selector<Scalar,int,Lower>::run); 175 176 init = true; 177 } 178 179 Scalar* x = reinterpret_cast<Scalar*>(px); 180 Scalar* y = reinterpret_cast<Scalar*>(py); 181 Scalar* ap = reinterpret_cast<Scalar*>(pap); 182 Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 183 184 int info = 0; 185 if(UPLO(*uplo)==INVALID) info = 1; 186 else if(*n<0) info = 2; 187 else if(*incx==0) info = 5; 188 else if(*incy==0) info = 7; 189 if(info) 190 return xerbla_(SCALAR_SUFFIX_UP"HPR2 ",&info,6); 191 192 if(alpha==Scalar(0)) 193 return 1; 194 195 Scalar* x_cpy = get_compact_vector(x, *n, *incx); 196 Scalar* y_cpy = get_compact_vector(y, *n, *incy); 197 198 int code = UPLO(*uplo); 199 if(code>=2 || func[code]==0) 200 return 0; 201 202 func[code](*n, ap, x_cpy, y_cpy, alpha); 203 204 if(x_cpy!=x) delete[] x_cpy; 205 if(y_cpy!=y) delete[] y_cpy; 206 207 return 1; 208 } 209 210 /** ZHER performs the hermitian rank 1 operation 211 * 212 * A := alpha*x*conjg( x' ) + A, 213 * 214 * where alpha is a real scalar, x is an n element vector and A is an 215 * n by n hermitian matrix. 216 */ 217 int EIGEN_BLAS_FUNC(her)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *pa, int *lda) 218 { 219 typedef void (*functype)(int, Scalar*, int, const Scalar*, const Scalar*, const Scalar&); 220 static functype func[2]; 221 222 static bool init = false; 223 if(!init) 224 { 225 for(int k=0; k<2; ++k) 226 func[k] = 0; 227 228 func[UP] = (selfadjoint_rank1_update<Scalar,int,ColMajor,Upper,false,Conj>::run); 229 func[LO] = (selfadjoint_rank1_update<Scalar,int,ColMajor,Lower,false,Conj>::run); 230 231 init = true; 232 } 233 234 Scalar* x = reinterpret_cast<Scalar*>(px); 235 Scalar* a = reinterpret_cast<Scalar*>(pa); 236 RealScalar alpha = *reinterpret_cast<RealScalar*>(palpha); 237 238 int info = 0; 239 if(UPLO(*uplo)==INVALID) info = 1; 240 else if(*n<0) info = 2; 241 else if(*incx==0) info = 5; 242 else if(*lda<std::max(1,*n)) info = 7; 243 if(info) 244 return xerbla_(SCALAR_SUFFIX_UP"HER ",&info,6); 245 246 if(alpha==RealScalar(0)) 247 return 1; 248 249 Scalar* x_cpy = get_compact_vector(x, *n, *incx); 250 251 int code = UPLO(*uplo); 252 if(code>=2 || func[code]==0) 253 return 0; 254 255 func[code](*n, a, *lda, x_cpy, x_cpy, alpha); 256 257 matrix(a,*n,*n,*lda).diagonal().imag().setZero(); 258 259 if(x_cpy!=x) delete[] x_cpy; 260 261 return 1; 262 } 263 264 /** ZHER2 performs the hermitian rank 2 operation 265 * 266 * A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A, 267 * 268 * where alpha is a scalar, x and y are n element vectors and A is an n 269 * by n hermitian matrix. 270 */ 271 int EIGEN_BLAS_FUNC(her2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pa, int *lda) 272 { 273 typedef void (*functype)(int, Scalar*, int, const Scalar*, const Scalar*, Scalar); 274 static functype func[2]; 275 276 static bool init = false; 277 if(!init) 278 { 279 for(int k=0; k<2; ++k) 280 func[k] = 0; 281 282 func[UP] = (internal::rank2_update_selector<Scalar,int,Upper>::run); 283 func[LO] = (internal::rank2_update_selector<Scalar,int,Lower>::run); 284 285 init = true; 286 } 287 288 Scalar* x = reinterpret_cast<Scalar*>(px); 289 Scalar* y = reinterpret_cast<Scalar*>(py); 290 Scalar* a = reinterpret_cast<Scalar*>(pa); 291 Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 292 293 int info = 0; 294 if(UPLO(*uplo)==INVALID) info = 1; 295 else if(*n<0) info = 2; 296 else if(*incx==0) info = 5; 297 else if(*incy==0) info = 7; 298 else if(*lda<std::max(1,*n)) info = 9; 299 if(info) 300 return xerbla_(SCALAR_SUFFIX_UP"HER2 ",&info,6); 301 302 if(alpha==Scalar(0)) 303 return 1; 304 305 Scalar* x_cpy = get_compact_vector(x, *n, *incx); 306 Scalar* y_cpy = get_compact_vector(y, *n, *incy); 307 308 int code = UPLO(*uplo); 309 if(code>=2 || func[code]==0) 310 return 0; 311 312 func[code](*n, a, *lda, x_cpy, y_cpy, alpha); 313 314 matrix(a,*n,*n,*lda).diagonal().imag().setZero(); 315 316 if(x_cpy!=x) delete[] x_cpy; 317 if(y_cpy!=y) delete[] y_cpy; 318 319 return 1; 320 } 321 322 /** ZGERU performs the rank 1 operation 323 * 324 * A := alpha*x*y' + A, 325 * 326 * where alpha is a scalar, x is an m element vector, y is an n element 327 * vector and A is an m by n matrix. 328 */ 329 int EIGEN_BLAS_FUNC(geru)(int *m, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pa, int *lda) 330 { 331 Scalar* x = reinterpret_cast<Scalar*>(px); 332 Scalar* y = reinterpret_cast<Scalar*>(py); 333 Scalar* a = reinterpret_cast<Scalar*>(pa); 334 Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 335 336 int info = 0; 337 if(*m<0) info = 1; 338 else if(*n<0) info = 2; 339 else if(*incx==0) info = 5; 340 else if(*incy==0) info = 7; 341 else if(*lda<std::max(1,*m)) info = 9; 342 if(info) 343 return xerbla_(SCALAR_SUFFIX_UP"GERU ",&info,6); 344 345 if(alpha==Scalar(0)) 346 return 1; 347 348 Scalar* x_cpy = get_compact_vector(x,*m,*incx); 349 Scalar* y_cpy = get_compact_vector(y,*n,*incy); 350 351 internal::general_rank1_update<Scalar,int,ColMajor,false,false>::run(*m, *n, a, *lda, x_cpy, y_cpy, alpha); 352 353 if(x_cpy!=x) delete[] x_cpy; 354 if(y_cpy!=y) delete[] y_cpy; 355 356 return 1; 357 } 358 359 /** ZGERC performs the rank 1 operation 360 * 361 * A := alpha*x*conjg( y' ) + A, 362 * 363 * where alpha is a scalar, x is an m element vector, y is an n element 364 * vector and A is an m by n matrix. 365 */ 366 int EIGEN_BLAS_FUNC(gerc)(int *m, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pa, int *lda) 367 { 368 Scalar* x = reinterpret_cast<Scalar*>(px); 369 Scalar* y = reinterpret_cast<Scalar*>(py); 370 Scalar* a = reinterpret_cast<Scalar*>(pa); 371 Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 372 373 int info = 0; 374 if(*m<0) info = 1; 375 else if(*n<0) info = 2; 376 else if(*incx==0) info = 5; 377 else if(*incy==0) info = 7; 378 else if(*lda<std::max(1,*m)) info = 9; 379 if(info) 380 return xerbla_(SCALAR_SUFFIX_UP"GERC ",&info,6); 381 382 if(alpha==Scalar(0)) 383 return 1; 384 385 Scalar* x_cpy = get_compact_vector(x,*m,*incx); 386 Scalar* y_cpy = get_compact_vector(y,*n,*incy); 387 388 internal::general_rank1_update<Scalar,int,ColMajor,false,Conj>::run(*m, *n, a, *lda, x_cpy, y_cpy, alpha); 389 390 if(x_cpy!=x) delete[] x_cpy; 391 if(y_cpy!=y) delete[] y_cpy; 392 393 return 1; 394 } 395
