ground/gcs/src/libs/eigen/blas/level2_impl.h

   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 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 }