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