[136] | 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 |
|
---|
| 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 |
|
---|