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 int EIGEN_BLAS_FUNC(gemv)(char *opa, int *m, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *incb, RealScalar *pbeta, RealScalar *pc, int *incc) 13 { 14 typedef void (*functype)(int, int, const Scalar *, int, const Scalar *, int , Scalar *, int, Scalar); 15 static functype func[4]; 16 17 static bool init = false; 18 if(!init) 19 { 20 for(int k=0; k<4; ++k) 21 func[k] = 0; 22 23 func[NOTR] = (internal::general_matrix_vector_product<int,Scalar,ColMajor,false,Scalar,false>::run); 24 func[TR ] = (internal::general_matrix_vector_product<int,Scalar,RowMajor,false,Scalar,false>::run); 25 func[ADJ ] = (internal::general_matrix_vector_product<int,Scalar,RowMajor,Conj, Scalar,false>::run); 26 27 init = true; 28 } 29 30 Scalar* a = reinterpret_cast<Scalar*>(pa); 31 Scalar* b = reinterpret_cast<Scalar*>(pb); 32 Scalar* c = reinterpret_cast<Scalar*>(pc); 33 Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 34 Scalar beta = *reinterpret_cast<Scalar*>(pbeta); 35 36 // check arguments 37 int info = 0; 38 if(OP(*opa)==INVALID) info = 1; 39 else if(*m<0) info = 2; 40 else if(*n<0) info = 3; 41 else if(*lda<std::max(1,*m)) info = 6; 42 else if(*incb==0) info = 8; 43 else if(*incc==0) info = 11; 44 if(info) 45 return xerbla_(SCALAR_SUFFIX_UP"GEMV ",&info,6); 46 47 if(*m==0 || *n==0 || (alpha==Scalar(0) && beta==Scalar(1))) 48 return 0; 49 50 int actual_m = *m; 51 int actual_n = *n; 52 int code = OP(*opa); 53 if(code!=NOTR) 54 std::swap(actual_m,actual_n); 55 56 Scalar* actual_b = get_compact_vector(b,actual_n,*incb); 57 Scalar* actual_c = get_compact_vector(c,actual_m,*incc); 58 59 if(beta!=Scalar(1)) 60 { 61 if(beta==Scalar(0)) vector(actual_c, actual_m).setZero(); 62 else vector(actual_c, actual_m) *= beta; 63 } 64 65 if(code>=4 || func[code]==0) 66 return 0; 67 68 func[code](actual_m, actual_n, a, *lda, actual_b, 1, actual_c, 1, alpha); 69 70 if(actual_b!=b) delete[] actual_b; 71 if(actual_c!=c) delete[] copy_back(actual_c,c,actual_m,*incc); 72 73 return 1; 74 } 75 76 int EIGEN_BLAS_FUNC(trsv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pa, int *lda, RealScalar *pb, int *incb) 77 { 78 typedef void (*functype)(int, const Scalar *, int, Scalar *); 79 static functype func[16]; 80 81 static bool init = false; 82 if(!init) 83 { 84 for(int k=0; k<16; ++k) 85 func[k] = 0; 86 87 func[NOTR | (UP << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0, false,ColMajor>::run); 88 func[TR | (UP << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0, false,RowMajor>::run); 89 func[ADJ | (UP << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0, Conj, RowMajor>::run); 90 91 func[NOTR | (LO << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0, false,ColMajor>::run); 92 func[TR | (LO << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0, false,RowMajor>::run); 93 func[ADJ | (LO << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0, Conj, RowMajor>::run); 94 95 func[NOTR | (UP << 2) | (UNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,ColMajor>::run); 96 func[TR | (UP << 2) | (UNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,RowMajor>::run); 97 func[ADJ | (UP << 2) | (UNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,Conj, RowMajor>::run); 98 99 func[NOTR | (LO << 2) | (UNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,ColMajor>::run); 100 func[TR | (LO << 2) | (UNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,RowMajor>::run); 101 func[ADJ | (LO << 2) | (UNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,Conj, RowMajor>::run); 102 103 init = true; 104 } 105 106 Scalar* a = reinterpret_cast<Scalar*>(pa); 107 Scalar* b = reinterpret_cast<Scalar*>(pb); 108 109 int info = 0; 110 if(UPLO(*uplo)==INVALID) info = 1; 111 else if(OP(*opa)==INVALID) info = 2; 112 else if(DIAG(*diag)==INVALID) info = 3; 113 else if(*n<0) info = 4; 114 else if(*lda<std::max(1,*n)) info = 6; 115 else if(*incb==0) info = 8; 116 if(info) 117 return xerbla_(SCALAR_SUFFIX_UP"TRSV ",&info,6); 118 119 Scalar* actual_b = get_compact_vector(b,*n,*incb); 120 121 int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3); 122 func[code](*n, a, *lda, actual_b); 123 124 if(actual_b!=b) delete[] copy_back(actual_b,b,*n,*incb); 125 126 return 0; 127 } 128 129 130 131 int EIGEN_BLAS_FUNC(trmv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pa, int *lda, RealScalar *pb, int *incb) 132 { 133 typedef void (*functype)(int, int, const Scalar *, int, const Scalar *, int, Scalar *, int, const Scalar&); 134 static functype func[16]; 135 136 static bool init = false; 137 if(!init) 138 { 139 for(int k=0; k<16; ++k) 140 func[k] = 0; 141 142 func[NOTR | (UP << 2) | (NUNIT << 3)] = (internal::triangular_matrix_vector_product<int,Upper|0, Scalar,false,Scalar,false,ColMajor>::run); 143 func[TR | (UP << 2) | (NUNIT << 3)] = (internal::triangular_matrix_vector_product<int,Lower|0, Scalar,false,Scalar,false,RowMajor>::run); 144 func[ADJ | (UP << 2) | (NUNIT << 3)] = (internal::triangular_matrix_vector_product<int,Lower|0, Scalar,Conj, Scalar,false,RowMajor>::run); 145 146 func[NOTR | (LO << 2) | (NUNIT << 3)] = (internal::triangular_matrix_vector_product<int,Lower|0, Scalar,false,Scalar,false,ColMajor>::run); 147 func[TR | (LO << 2) | (NUNIT << 3)] = (internal::triangular_matrix_vector_product<int,Upper|0, Scalar,false,Scalar,false,RowMajor>::run); 148 func[ADJ | (LO << 2) | (NUNIT << 3)] = (internal::triangular_matrix_vector_product<int,Upper|0, Scalar,Conj, Scalar,false,RowMajor>::run); 149 150 func[NOTR | (UP << 2) | (UNIT << 3)] = (internal::triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run); 151 func[TR | (UP << 2) | (UNIT << 3)] = (internal::triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run); 152 func[ADJ | (UP << 2) | (UNIT << 3)] = (internal::triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run); 153 154 func[NOTR | (LO << 2) | (UNIT << 3)] = (internal::triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run); 155 func[TR | (LO << 2) | (UNIT << 3)] = (internal::triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run); 156 func[ADJ | (LO << 2) | (UNIT << 3)] = (internal::triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run); 157 158 init = true; 159 } 160 161 Scalar* a = reinterpret_cast<Scalar*>(pa); 162 Scalar* b = reinterpret_cast<Scalar*>(pb); 163 164 int info = 0; 165 if(UPLO(*uplo)==INVALID) info = 1; 166 else if(OP(*opa)==INVALID) info = 2; 167 else if(DIAG(*diag)==INVALID) info = 3; 168 else if(*n<0) info = 4; 169 else if(*lda<std::max(1,*n)) info = 6; 170 else if(*incb==0) info = 8; 171 if(info) 172 return xerbla_(SCALAR_SUFFIX_UP"TRMV ",&info,6); 173 174 if(*n==0) 175 return 1; 176 177 Scalar* actual_b = get_compact_vector(b,*n,*incb); 178 Matrix<Scalar,Dynamic,1> res(*n); 179 res.setZero(); 180 181 int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3); 182 if(code>=16 || func[code]==0) 183 return 0; 184 185 func[code](*n, *n, a, *lda, actual_b, 1, res.data(), 1, Scalar(1)); 186 187 copy_back(res.data(),b,*n,*incb); 188 if(actual_b!=b) delete[] actual_b; 189 190 return 1; 191 } 192 193 /** GBMV performs one of the matrix-vector operations 194 * 195 * y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, 196 * 197 * where alpha and beta are scalars, x and y are vectors and A is an 198 * m by n band matrix, with kl sub-diagonals and ku super-diagonals. 199 */ 200 int EIGEN_BLAS_FUNC(gbmv)(char *trans, int *m, int *n, int *kl, int *ku, RealScalar *palpha, RealScalar *pa, int *lda, 201 RealScalar *px, int *incx, RealScalar *pbeta, RealScalar *py, int *incy) 202 { 203 Scalar* a = reinterpret_cast<Scalar*>(pa); 204 Scalar* x = reinterpret_cast<Scalar*>(px); 205 Scalar* y = reinterpret_cast<Scalar*>(py); 206 Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 207 Scalar beta = *reinterpret_cast<Scalar*>(pbeta); 208 int coeff_rows = *kl+*ku+1; 209 210 int info = 0; 211 if(OP(*trans)==INVALID) info = 1; 212 else if(*m<0) info = 2; 213 else if(*n<0) info = 3; 214 else if(*kl<0) info = 4; 215 else if(*ku<0) info = 5; 216 else if(*lda<coeff_rows) info = 8; 217 else if(*incx==0) info = 10; 218 else if(*incy==0) info = 13; 219 if(info) 220 return xerbla_(SCALAR_SUFFIX_UP"GBMV ",&info,6); 221 222 if(*m==0 || *n==0 || (alpha==Scalar(0) && beta==Scalar(1))) 223 return 0; 224 225 int actual_m = *m; 226 int actual_n = *n; 227 if(OP(*trans)!=NOTR) 228 std::swap(actual_m,actual_n); 229 230 Scalar* actual_x = get_compact_vector(x,actual_n,*incx); 231 Scalar* actual_y = get_compact_vector(y,actual_m,*incy); 232 233 if(beta!=Scalar(1)) 234 { 235 if(beta==Scalar(0)) vector(actual_y, actual_m).setZero(); 236 else vector(actual_y, actual_m) *= beta; 237 } 238 239 MatrixType mat_coeffs(a,coeff_rows,*n,*lda); 240 241 int nb = std::min(*n,(*m)+(*ku)); 242 for(int j=0; j<nb; ++j) 243 { 244 int start = std::max(0,j - *ku); 245 int end = std::min((*m)-1,j + *kl); 246 int len = end - start + 1; 247 int offset = (*ku) - j + start; 248 if(OP(*trans)==NOTR) 249 vector(actual_y+start,len) += (alpha*actual_x[j]) * mat_coeffs.col(j).segment(offset,len); 250 else if(OP(*trans)==TR) 251 actual_y[j] += alpha * ( mat_coeffs.col(j).segment(offset,len).transpose() * vector(actual_x+start,len) ).value(); 252 else 253 actual_y[j] += alpha * ( mat_coeffs.col(j).segment(offset,len).adjoint() * vector(actual_x+start,len) ).value(); 254 } 255 256 if(actual_x!=x) delete[] actual_x; 257 if(actual_y!=y) delete[] copy_back(actual_y,y,actual_m,*incy); 258 259 return 0; 260 } 261 262 #if 0 263 /** TBMV performs one of the matrix-vector operations 264 * 265 * x := A*x, or x := A'*x, 266 * 267 * where x is an n element vector and A is an n by n unit, or non-unit, 268 * upper or lower triangular band matrix, with ( k + 1 ) diagonals. 269 */ 270 int EIGEN_BLAS_FUNC(tbmv)(char *uplo, char *opa, char *diag, int *n, int *k, RealScalar *pa, int *lda, RealScalar *px, int *incx) 271 { 272 Scalar* a = reinterpret_cast<Scalar*>(pa); 273 Scalar* x = reinterpret_cast<Scalar*>(px); 274 int coeff_rows = *k + 1; 275 276 int info = 0; 277 if(UPLO(*uplo)==INVALID) info = 1; 278 else if(OP(*opa)==INVALID) info = 2; 279 else if(DIAG(*diag)==INVALID) info = 3; 280 else if(*n<0) info = 4; 281 else if(*k<0) info = 5; 282 else if(*lda<coeff_rows) info = 7; 283 else if(*incx==0) info = 9; 284 if(info) 285 return xerbla_(SCALAR_SUFFIX_UP"TBMV ",&info,6); 286 287 if(*n==0) 288 return 0; 289 290 int actual_n = *n; 291 292 Scalar* actual_x = get_compact_vector(x,actual_n,*incx); 293 294 MatrixType mat_coeffs(a,coeff_rows,*n,*lda); 295 296 int ku = UPLO(*uplo)==UPPER ? *k : 0; 297 int kl = UPLO(*uplo)==LOWER ? *k : 0; 298 299 for(int j=0; j<*n; ++j) 300 { 301 int start = std::max(0,j - ku); 302 int end = std::min((*m)-1,j + kl); 303 int len = end - start + 1; 304 int offset = (ku) - j + start; 305 306 if(OP(*trans)==NOTR) 307 vector(actual_y+start,len) += (alpha*actual_x[j]) * mat_coeffs.col(j).segment(offset,len); 308 else if(OP(*trans)==TR) 309 actual_y[j] += alpha * ( mat_coeffs.col(j).segment(offset,len).transpose() * vector(actual_x+start,len) ).value(); 310 else 311 actual_y[j] += alpha * ( mat_coeffs.col(j).segment(offset,len).adjoint() * vector(actual_x+start,len) ).value(); 312 } 313 314 if(actual_x!=x) delete[] actual_x; 315 if(actual_y!=y) delete[] copy_back(actual_y,y,actual_m,*incy); 316 317 return 0; 318 } 319 #endif 320 321 /** DTBSV solves one of the systems of equations 322 * 323 * A*x = b, or A'*x = b, 324 * 325 * where b and x are n element vectors and A is an n by n unit, or 326 * non-unit, upper or lower triangular band matrix, with ( k + 1 ) 327 * diagonals. 328 * 329 * No test for singularity or near-singularity is included in this 330 * routine. Such tests must be performed before calling this routine. 331 */ 332 int EIGEN_BLAS_FUNC(tbsv)(char *uplo, char *op, char *diag, int *n, int *k, RealScalar *pa, int *lda, RealScalar *px, int *incx) 333 { 334 typedef void (*functype)(int, int, const Scalar *, int, Scalar *); 335 static functype func[16]; 336 337 static bool init = false; 338 if(!init) 339 { 340 for(int k=0; k<16; ++k) 341 func[k] = 0; 342 343 func[NOTR | (UP << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector<int,Upper|0, Scalar,false,Scalar,ColMajor>::run); 344 func[TR | (UP << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector<int,Lower|0, Scalar,false,Scalar,RowMajor>::run); 345 func[ADJ | (UP << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector<int,Lower|0, Scalar,Conj, Scalar,RowMajor>::run); 346 347 func[NOTR | (LO << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector<int,Lower|0, Scalar,false,Scalar,ColMajor>::run); 348 func[TR | (LO << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector<int,Upper|0, Scalar,false,Scalar,RowMajor>::run); 349 func[ADJ | (LO << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector<int,Upper|0, Scalar,Conj, Scalar,RowMajor>::run); 350 351 func[NOTR | (UP << 2) | (UNIT << 3)] = (internal::band_solve_triangular_selector<int,Upper|UnitDiag,Scalar,false,Scalar,ColMajor>::run); 352 func[TR | (UP << 2) | (UNIT << 3)] = (internal::band_solve_triangular_selector<int,Lower|UnitDiag,Scalar,false,Scalar,RowMajor>::run); 353 func[ADJ | (UP << 2) | (UNIT << 3)] = (internal::band_solve_triangular_selector<int,Lower|UnitDiag,Scalar,Conj, Scalar,RowMajor>::run); 354 355 func[NOTR | (LO << 2) | (UNIT << 3)] = (internal::band_solve_triangular_selector<int,Lower|UnitDiag,Scalar,false,Scalar,ColMajor>::run); 356 func[TR | (LO << 2) | (UNIT << 3)] = (internal::band_solve_triangular_selector<int,Upper|UnitDiag,Scalar,false,Scalar,RowMajor>::run); 357 func[ADJ | (LO << 2) | (UNIT << 3)] = (internal::band_solve_triangular_selector<int,Upper|UnitDiag,Scalar,Conj, Scalar,RowMajor>::run); 358 359 init = true; 360 } 361 362 Scalar* a = reinterpret_cast<Scalar*>(pa); 363 Scalar* x = reinterpret_cast<Scalar*>(px); 364 int coeff_rows = *k+1; 365 366 int info = 0; 367 if(UPLO(*uplo)==INVALID) info = 1; 368 else if(OP(*op)==INVALID) info = 2; 369 else if(DIAG(*diag)==INVALID) info = 3; 370 else if(*n<0) info = 4; 371 else if(*k<0) info = 5; 372 else if(*lda<coeff_rows) info = 7; 373 else if(*incx==0) info = 9; 374 if(info) 375 return xerbla_(SCALAR_SUFFIX_UP"TBSV ",&info,6); 376 377 if(*n==0 || (*k==0 && DIAG(*diag)==UNIT)) 378 return 0; 379 380 int actual_n = *n; 381 382 Scalar* actual_x = get_compact_vector(x,actual_n,*incx); 383 384 int code = OP(*op) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3); 385 if(code>=16 || func[code]==0) 386 return 0; 387 388 func[code](*n, *k, a, *lda, actual_x); 389 390 if(actual_x!=x) delete[] copy_back(actual_x,x,actual_n,*incx); 391 392 return 0; 393 } 394 395 /** DTPMV performs one of the matrix-vector operations 396 * 397 * x := A*x, or x := A'*x, 398 * 399 * where x is an n element vector and A is an n by n unit, or non-unit, 400 * upper or lower triangular matrix, supplied in packed form. 401 */ 402 int EIGEN_BLAS_FUNC(tpmv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pap, RealScalar *px, int *incx) 403 { 404 typedef void (*functype)(int, const Scalar*, const Scalar*, Scalar*, Scalar); 405 static functype func[16]; 406 407 static bool init = false; 408 if(!init) 409 { 410 for(int k=0; k<16; ++k) 411 func[k] = 0; 412 413 func[NOTR | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|0, Scalar,false,Scalar,false,ColMajor>::run); 414 func[TR | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|0, Scalar,false,Scalar,false,RowMajor>::run); 415 func[ADJ | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|0, Scalar,Conj, Scalar,false,RowMajor>::run); 416 417 func[NOTR | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|0, Scalar,false,Scalar,false,ColMajor>::run); 418 func[TR | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|0, Scalar,false,Scalar,false,RowMajor>::run); 419 func[ADJ | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|0, Scalar,Conj, Scalar,false,RowMajor>::run); 420 421 func[NOTR | (UP << 2) | (UNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run); 422 func[TR | (UP << 2) | (UNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run); 423 func[ADJ | (UP << 2) | (UNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run); 424 425 func[NOTR | (LO << 2) | (UNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Lower|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run); 426 func[TR | (LO << 2) | (UNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run); 427 func[ADJ | (LO << 2) | (UNIT << 3)] = (internal::packed_triangular_matrix_vector_product<int,Upper|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run); 428 429 init = true; 430 } 431 432 Scalar* ap = reinterpret_cast<Scalar*>(pap); 433 Scalar* x = reinterpret_cast<Scalar*>(px); 434 435 int info = 0; 436 if(UPLO(*uplo)==INVALID) info = 1; 437 else if(OP(*opa)==INVALID) info = 2; 438 else if(DIAG(*diag)==INVALID) info = 3; 439 else if(*n<0) info = 4; 440 else if(*incx==0) info = 7; 441 if(info) 442 return xerbla_(SCALAR_SUFFIX_UP"TPMV ",&info,6); 443 444 if(*n==0) 445 return 1; 446 447 Scalar* actual_x = get_compact_vector(x,*n,*incx); 448 Matrix<Scalar,Dynamic,1> res(*n); 449 res.setZero(); 450 451 int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3); 452 if(code>=16 || func[code]==0) 453 return 0; 454 455 func[code](*n, ap, actual_x, res.data(), Scalar(1)); 456 457 copy_back(res.data(),x,*n,*incx); 458 if(actual_x!=x) delete[] actual_x; 459 460 return 1; 461 } 462 463 /** DTPSV solves one of the systems of equations 464 * 465 * A*x = b, or A'*x = b, 466 * 467 * where b and x are n element vectors and A is an n by n unit, or 468 * non-unit, upper or lower triangular matrix, supplied in packed form. 469 * 470 * No test for singularity or near-singularity is included in this 471 * routine. Such tests must be performed before calling this routine. 472 */ 473 int EIGEN_BLAS_FUNC(tpsv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pap, RealScalar *px, int *incx) 474 { 475 typedef void (*functype)(int, const Scalar*, Scalar*); 476 static functype func[16]; 477 478 static bool init = false; 479 if(!init) 480 { 481 for(int k=0; k<16; ++k) 482 func[k] = 0; 483 484 func[NOTR | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0, false,ColMajor>::run); 485 func[TR | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0, false,RowMajor>::run); 486 func[ADJ | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0, Conj, RowMajor>::run); 487 488 func[NOTR | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0, false,ColMajor>::run); 489 func[TR | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0, false,RowMajor>::run); 490 func[ADJ | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0, Conj, RowMajor>::run); 491 492 func[NOTR | (UP << 2) | (UNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,ColMajor>::run); 493 func[TR | (UP << 2) | (UNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,RowMajor>::run); 494 func[ADJ | (UP << 2) | (UNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,Conj, RowMajor>::run); 495 496 func[NOTR | (LO << 2) | (UNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,ColMajor>::run); 497 func[TR | (LO << 2) | (UNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,RowMajor>::run); 498 func[ADJ | (LO << 2) | (UNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,Conj, RowMajor>::run); 499 500 init = true; 501 } 502 503 Scalar* ap = reinterpret_cast<Scalar*>(pap); 504 Scalar* x = reinterpret_cast<Scalar*>(px); 505 506 int info = 0; 507 if(UPLO(*uplo)==INVALID) info = 1; 508 else if(OP(*opa)==INVALID) info = 2; 509 else if(DIAG(*diag)==INVALID) info = 3; 510 else if(*n<0) info = 4; 511 else if(*incx==0) info = 7; 512 if(info) 513 return xerbla_(SCALAR_SUFFIX_UP"TPSV ",&info,6); 514 515 Scalar* actual_x = get_compact_vector(x,*n,*incx); 516 517 int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3); 518 func[code](*n, ap, actual_x); 519 520 if(actual_x!=x) delete[] copy_back(actual_x,x,*n,*incx); 521 522 return 1; 523 } 524 525
