level2_impl.h 25 KB

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