level2_cplx_impl.h 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360
  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. /** ZHEMV performs the matrix-vector operation
  11. *
  12. * y := alpha*A*x + beta*y,
  13. *
  14. * where alpha and beta are scalars, x and y are n element vectors and
  15. * A is an n by n hermitian matrix.
  16. */
  17. int EIGEN_BLAS_FUNC(hemv)(const char *uplo, const int *n, const RealScalar *palpha, const RealScalar *pa, const int *lda,
  18. const RealScalar *px, const int *incx, const RealScalar *pbeta, RealScalar *py, const int *incy)
  19. {
  20. typedef void (*functype)(int, const Scalar*, int, const Scalar*, Scalar*, Scalar);
  21. static const functype func[2] = {
  22. // array index: UP
  23. (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Upper,false,false>::run),
  24. // array index: LO
  25. (internal::selfadjoint_matrix_vector_product<Scalar,int,ColMajor,Lower,false,false>::run),
  26. };
  27. const Scalar* a = reinterpret_cast<const Scalar*>(pa);
  28. const Scalar* x = reinterpret_cast<const Scalar*>(px);
  29. Scalar* y = reinterpret_cast<Scalar*>(py);
  30. Scalar alpha = *reinterpret_cast<const Scalar*>(palpha);
  31. Scalar beta = *reinterpret_cast<const Scalar*>(pbeta);
  32. // check arguments
  33. int info = 0;
  34. if(UPLO(*uplo)==INVALID) info = 1;
  35. else if(*n<0) info = 2;
  36. else if(*lda<std::max(1,*n)) info = 5;
  37. else if(*incx==0) info = 7;
  38. else if(*incy==0) info = 10;
  39. if(info)
  40. return xerbla_(SCALAR_SUFFIX_UP"HEMV ",&info,6);
  41. if(*n==0)
  42. return 1;
  43. const Scalar* actual_x = get_compact_vector(x,*n,*incx);
  44. Scalar* actual_y = get_compact_vector(y,*n,*incy);
  45. if(beta!=Scalar(1))
  46. {
  47. if(beta==Scalar(0)) make_vector(actual_y, *n).setZero();
  48. else make_vector(actual_y, *n) *= beta;
  49. }
  50. if(alpha!=Scalar(0))
  51. {
  52. int code = UPLO(*uplo);
  53. if(code>=2 || func[code]==0)
  54. return 0;
  55. func[code](*n, a, *lda, actual_x, actual_y, alpha);
  56. }
  57. if(actual_x!=x) delete[] actual_x;
  58. if(actual_y!=y) delete[] copy_back(actual_y,y,*n,*incy);
  59. return 1;
  60. }
  61. /** ZHBMV performs the matrix-vector operation
  62. *
  63. * y := alpha*A*x + beta*y,
  64. *
  65. * where alpha and beta are scalars, x and y are n element vectors and
  66. * A is an n by n hermitian band matrix, with k super-diagonals.
  67. */
  68. // int EIGEN_BLAS_FUNC(hbmv)(char *uplo, int *n, int *k, RealScalar *alpha, RealScalar *a, int *lda,
  69. // RealScalar *x, int *incx, RealScalar *beta, RealScalar *y, int *incy)
  70. // {
  71. // return 1;
  72. // }
  73. /** ZHPMV performs the matrix-vector operation
  74. *
  75. * y := alpha*A*x + beta*y,
  76. *
  77. * where alpha and beta are scalars, x and y are n element vectors and
  78. * A is an n by n hermitian matrix, supplied in packed form.
  79. */
  80. // int EIGEN_BLAS_FUNC(hpmv)(char *uplo, int *n, RealScalar *alpha, RealScalar *ap, RealScalar *x, int *incx, RealScalar *beta, RealScalar *y, int *incy)
  81. // {
  82. // return 1;
  83. // }
  84. /** ZHPR performs the hermitian rank 1 operation
  85. *
  86. * A := alpha*x*conjg( x' ) + A,
  87. *
  88. * where alpha is a real scalar, x is an n element vector and A is an
  89. * n by n hermitian matrix, supplied in packed form.
  90. */
  91. int EIGEN_BLAS_FUNC(hpr)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *pap)
  92. {
  93. typedef void (*functype)(int, Scalar*, const Scalar*, RealScalar);
  94. static const functype func[2] = {
  95. // array index: UP
  96. (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Upper,false,Conj>::run),
  97. // array index: LO
  98. (internal::selfadjoint_packed_rank1_update<Scalar,int,ColMajor,Lower,false,Conj>::run),
  99. };
  100. Scalar* x = reinterpret_cast<Scalar*>(px);
  101. Scalar* ap = reinterpret_cast<Scalar*>(pap);
  102. RealScalar alpha = *palpha;
  103. int info = 0;
  104. if(UPLO(*uplo)==INVALID) info = 1;
  105. else if(*n<0) info = 2;
  106. else if(*incx==0) info = 5;
  107. if(info)
  108. return xerbla_(SCALAR_SUFFIX_UP"HPR ",&info,6);
  109. if(alpha==Scalar(0))
  110. return 1;
  111. Scalar* x_cpy = get_compact_vector(x, *n, *incx);
  112. int code = UPLO(*uplo);
  113. if(code>=2 || func[code]==0)
  114. return 0;
  115. func[code](*n, ap, x_cpy, alpha);
  116. if(x_cpy!=x) delete[] x_cpy;
  117. return 1;
  118. }
  119. /** ZHPR2 performs the hermitian rank 2 operation
  120. *
  121. * A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A,
  122. *
  123. * where alpha is a scalar, x and y are n element vectors and A is an
  124. * n by n hermitian matrix, supplied in packed form.
  125. */
  126. int EIGEN_BLAS_FUNC(hpr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pap)
  127. {
  128. typedef void (*functype)(int, Scalar*, const Scalar*, const Scalar*, Scalar);
  129. static const functype func[2] = {
  130. // array index: UP
  131. (internal::packed_rank2_update_selector<Scalar,int,Upper>::run),
  132. // array index: LO
  133. (internal::packed_rank2_update_selector<Scalar,int,Lower>::run),
  134. };
  135. Scalar* x = reinterpret_cast<Scalar*>(px);
  136. Scalar* y = reinterpret_cast<Scalar*>(py);
  137. Scalar* ap = reinterpret_cast<Scalar*>(pap);
  138. Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
  139. int info = 0;
  140. if(UPLO(*uplo)==INVALID) info = 1;
  141. else if(*n<0) info = 2;
  142. else if(*incx==0) info = 5;
  143. else if(*incy==0) info = 7;
  144. if(info)
  145. return xerbla_(SCALAR_SUFFIX_UP"HPR2 ",&info,6);
  146. if(alpha==Scalar(0))
  147. return 1;
  148. Scalar* x_cpy = get_compact_vector(x, *n, *incx);
  149. Scalar* y_cpy = get_compact_vector(y, *n, *incy);
  150. int code = UPLO(*uplo);
  151. if(code>=2 || func[code]==0)
  152. return 0;
  153. func[code](*n, ap, x_cpy, y_cpy, alpha);
  154. if(x_cpy!=x) delete[] x_cpy;
  155. if(y_cpy!=y) delete[] y_cpy;
  156. return 1;
  157. }
  158. /** ZHER performs the hermitian rank 1 operation
  159. *
  160. * A := alpha*x*conjg( x' ) + A,
  161. *
  162. * where alpha is a real scalar, x is an n element vector and A is an
  163. * n by n hermitian matrix.
  164. */
  165. int EIGEN_BLAS_FUNC(her)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *pa, int *lda)
  166. {
  167. typedef void (*functype)(int, Scalar*, int, const Scalar*, const Scalar*, const Scalar&);
  168. static const functype func[2] = {
  169. // array index: UP
  170. (selfadjoint_rank1_update<Scalar,int,ColMajor,Upper,false,Conj>::run),
  171. // array index: LO
  172. (selfadjoint_rank1_update<Scalar,int,ColMajor,Lower,false,Conj>::run),
  173. };
  174. Scalar* x = reinterpret_cast<Scalar*>(px);
  175. Scalar* a = reinterpret_cast<Scalar*>(pa);
  176. RealScalar alpha = *reinterpret_cast<RealScalar*>(palpha);
  177. int info = 0;
  178. if(UPLO(*uplo)==INVALID) info = 1;
  179. else if(*n<0) info = 2;
  180. else if(*incx==0) info = 5;
  181. else if(*lda<std::max(1,*n)) info = 7;
  182. if(info)
  183. return xerbla_(SCALAR_SUFFIX_UP"HER ",&info,6);
  184. if(alpha==RealScalar(0))
  185. return 1;
  186. Scalar* x_cpy = get_compact_vector(x, *n, *incx);
  187. int code = UPLO(*uplo);
  188. if(code>=2 || func[code]==0)
  189. return 0;
  190. func[code](*n, a, *lda, x_cpy, x_cpy, alpha);
  191. matrix(a,*n,*n,*lda).diagonal().imag().setZero();
  192. if(x_cpy!=x) delete[] x_cpy;
  193. return 1;
  194. }
  195. /** ZHER2 performs the hermitian rank 2 operation
  196. *
  197. * A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A,
  198. *
  199. * where alpha is a scalar, x and y are n element vectors and A is an n
  200. * by n hermitian matrix.
  201. */
  202. int EIGEN_BLAS_FUNC(her2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pa, int *lda)
  203. {
  204. typedef void (*functype)(int, Scalar*, int, const Scalar*, const Scalar*, Scalar);
  205. static const functype func[2] = {
  206. // array index: UP
  207. (internal::rank2_update_selector<Scalar,int,Upper>::run),
  208. // array index: LO
  209. (internal::rank2_update_selector<Scalar,int,Lower>::run),
  210. };
  211. Scalar* x = reinterpret_cast<Scalar*>(px);
  212. Scalar* y = reinterpret_cast<Scalar*>(py);
  213. Scalar* a = reinterpret_cast<Scalar*>(pa);
  214. Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
  215. int info = 0;
  216. if(UPLO(*uplo)==INVALID) info = 1;
  217. else if(*n<0) info = 2;
  218. else if(*incx==0) info = 5;
  219. else if(*incy==0) info = 7;
  220. else if(*lda<std::max(1,*n)) info = 9;
  221. if(info)
  222. return xerbla_(SCALAR_SUFFIX_UP"HER2 ",&info,6);
  223. if(alpha==Scalar(0))
  224. return 1;
  225. Scalar* x_cpy = get_compact_vector(x, *n, *incx);
  226. Scalar* y_cpy = get_compact_vector(y, *n, *incy);
  227. int code = UPLO(*uplo);
  228. if(code>=2 || func[code]==0)
  229. return 0;
  230. func[code](*n, a, *lda, x_cpy, y_cpy, alpha);
  231. matrix(a,*n,*n,*lda).diagonal().imag().setZero();
  232. if(x_cpy!=x) delete[] x_cpy;
  233. if(y_cpy!=y) delete[] y_cpy;
  234. return 1;
  235. }
  236. /** ZGERU performs the rank 1 operation
  237. *
  238. * A := alpha*x*y' + A,
  239. *
  240. * where alpha is a scalar, x is an m element vector, y is an n element
  241. * vector and A is an m by n matrix.
  242. */
  243. int EIGEN_BLAS_FUNC(geru)(int *m, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pa, int *lda)
  244. {
  245. Scalar* x = reinterpret_cast<Scalar*>(px);
  246. Scalar* y = reinterpret_cast<Scalar*>(py);
  247. Scalar* a = reinterpret_cast<Scalar*>(pa);
  248. Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
  249. int info = 0;
  250. if(*m<0) info = 1;
  251. else if(*n<0) info = 2;
  252. else if(*incx==0) info = 5;
  253. else if(*incy==0) info = 7;
  254. else if(*lda<std::max(1,*m)) info = 9;
  255. if(info)
  256. return xerbla_(SCALAR_SUFFIX_UP"GERU ",&info,6);
  257. if(alpha==Scalar(0))
  258. return 1;
  259. Scalar* x_cpy = get_compact_vector(x,*m,*incx);
  260. Scalar* y_cpy = get_compact_vector(y,*n,*incy);
  261. internal::general_rank1_update<Scalar,int,ColMajor,false,false>::run(*m, *n, a, *lda, x_cpy, y_cpy, alpha);
  262. if(x_cpy!=x) delete[] x_cpy;
  263. if(y_cpy!=y) delete[] y_cpy;
  264. return 1;
  265. }
  266. /** ZGERC performs the rank 1 operation
  267. *
  268. * A := alpha*x*conjg( y' ) + A,
  269. *
  270. * where alpha is a scalar, x is an m element vector, y is an n element
  271. * vector and A is an m by n matrix.
  272. */
  273. int EIGEN_BLAS_FUNC(gerc)(int *m, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pa, int *lda)
  274. {
  275. Scalar* x = reinterpret_cast<Scalar*>(px);
  276. Scalar* y = reinterpret_cast<Scalar*>(py);
  277. Scalar* a = reinterpret_cast<Scalar*>(pa);
  278. Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
  279. int info = 0;
  280. if(*m<0) info = 1;
  281. else if(*n<0) info = 2;
  282. else if(*incx==0) info = 5;
  283. else if(*incy==0) info = 7;
  284. else if(*lda<std::max(1,*m)) info = 9;
  285. if(info)
  286. return xerbla_(SCALAR_SUFFIX_UP"GERC ",&info,6);
  287. if(alpha==Scalar(0))
  288. return 1;
  289. Scalar* x_cpy = get_compact_vector(x,*m,*incx);
  290. Scalar* y_cpy = get_compact_vector(y,*n,*incy);
  291. internal::general_rank1_update<Scalar,int,ColMajor,false,Conj>::run(*m, *n, a, *lda, x_cpy, y_cpy, alpha);
  292. if(x_cpy!=x) delete[] x_cpy;
  293. if(y_cpy!=y) delete[] y_cpy;
  294. return 1;
  295. }