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