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