/foreignlibs/boost/numeric/bindings/lapack/driver/ggevx.hpp

https://github.com/TRIQS/triqs_0.x · C++ Header · 755 lines · 581 code · 35 blank · 139 comment · 83 complexity · ac5df2909fc8e82c9003902ba5199e54 MD5 · raw file

  1. //
  2. // Copyright (c) 2002--2010
  3. // Toon Knapen, Karl Meerbergen, Kresimir Fresl,
  4. // Thomas Klimpel and Rutger ter Borg
  5. //
  6. // Distributed under the Boost Software License, Version 1.0.
  7. // (See accompanying file LICENSE_1_0.txt or copy at
  8. // http://www.boost.org/LICENSE_1_0.txt)
  9. //
  10. // THIS FILE IS AUTOMATICALLY GENERATED
  11. // PLEASE DO NOT EDIT!
  12. //
  13. #ifndef BOOST_NUMERIC_BINDINGS_LAPACK_DRIVER_GGEVX_HPP
  14. #define BOOST_NUMERIC_BINDINGS_LAPACK_DRIVER_GGEVX_HPP
  15. #include <boost/assert.hpp>
  16. #include <boost/numeric/bindings/begin.hpp>
  17. #include <boost/numeric/bindings/detail/array.hpp>
  18. #include <boost/numeric/bindings/is_column_major.hpp>
  19. #include <boost/numeric/bindings/is_complex.hpp>
  20. #include <boost/numeric/bindings/is_mutable.hpp>
  21. #include <boost/numeric/bindings/is_real.hpp>
  22. #include <boost/numeric/bindings/lapack/workspace.hpp>
  23. #include <boost/numeric/bindings/remove_imaginary.hpp>
  24. #include <boost/numeric/bindings/size.hpp>
  25. #include <boost/numeric/bindings/stride.hpp>
  26. #include <boost/numeric/bindings/traits/detail/utils.hpp>
  27. #include <boost/numeric/bindings/value_type.hpp>
  28. #include <boost/static_assert.hpp>
  29. #include <boost/type_traits/is_same.hpp>
  30. #include <boost/type_traits/remove_const.hpp>
  31. #include <boost/utility/enable_if.hpp>
  32. //
  33. // The LAPACK-backend for ggevx is the netlib-compatible backend.
  34. //
  35. #include <boost/numeric/bindings/lapack/detail/lapack.h>
  36. #include <boost/numeric/bindings/lapack/detail/lapack_option.hpp>
  37. namespace boost {
  38. namespace numeric {
  39. namespace bindings {
  40. namespace lapack {
  41. //
  42. // The detail namespace contains value-type-overloaded functions that
  43. // dispatch to the appropriate back-end LAPACK-routine.
  44. //
  45. namespace detail {
  46. //
  47. // Overloaded function for dispatching to
  48. // * netlib-compatible LAPACK backend (the default), and
  49. // * float value-type.
  50. //
  51. inline std::ptrdiff_t ggevx( const char balanc, const char jobvl,
  52. const char jobvr, const char sense, const fortran_int_t n, float* a,
  53. const fortran_int_t lda, float* b, const fortran_int_t ldb,
  54. float* alphar, float* alphai, float* beta, float* vl,
  55. const fortran_int_t ldvl, float* vr, const fortran_int_t ldvr,
  56. fortran_int_t& ilo, fortran_int_t& ihi, float* lscale, float* rscale,
  57. float& abnrm, float& bbnrm, float* rconde, float* rcondv, float* work,
  58. const fortran_int_t lwork, fortran_int_t* iwork,
  59. fortran_bool_t* bwork ) {
  60. fortran_int_t info(0);
  61. LAPACK_SGGEVX( &balanc, &jobvl, &jobvr, &sense, &n, a, &lda, b, &ldb,
  62. alphar, alphai, beta, vl, &ldvl, vr, &ldvr, &ilo, &ihi, lscale,
  63. rscale, &abnrm, &bbnrm, rconde, rcondv, work, &lwork, iwork,
  64. bwork, &info );
  65. return info;
  66. }
  67. //
  68. // Overloaded function for dispatching to
  69. // * netlib-compatible LAPACK backend (the default), and
  70. // * double value-type.
  71. //
  72. inline std::ptrdiff_t ggevx( const char balanc, const char jobvl,
  73. const char jobvr, const char sense, const fortran_int_t n, double* a,
  74. const fortran_int_t lda, double* b, const fortran_int_t ldb,
  75. double* alphar, double* alphai, double* beta, double* vl,
  76. const fortran_int_t ldvl, double* vr, const fortran_int_t ldvr,
  77. fortran_int_t& ilo, fortran_int_t& ihi, double* lscale,
  78. double* rscale, double& abnrm, double& bbnrm, double* rconde,
  79. double* rcondv, double* work, const fortran_int_t lwork,
  80. fortran_int_t* iwork, fortran_bool_t* bwork ) {
  81. fortran_int_t info(0);
  82. LAPACK_DGGEVX( &balanc, &jobvl, &jobvr, &sense, &n, a, &lda, b, &ldb,
  83. alphar, alphai, beta, vl, &ldvl, vr, &ldvr, &ilo, &ihi, lscale,
  84. rscale, &abnrm, &bbnrm, rconde, rcondv, work, &lwork, iwork,
  85. bwork, &info );
  86. return info;
  87. }
  88. //
  89. // Overloaded function for dispatching to
  90. // * netlib-compatible LAPACK backend (the default), and
  91. // * complex<float> value-type.
  92. //
  93. inline std::ptrdiff_t ggevx( const char balanc, const char jobvl,
  94. const char jobvr, const char sense, const fortran_int_t n,
  95. std::complex<float>* a, const fortran_int_t lda,
  96. std::complex<float>* b, const fortran_int_t ldb,
  97. std::complex<float>* alpha, std::complex<float>* beta,
  98. std::complex<float>* vl, const fortran_int_t ldvl,
  99. std::complex<float>* vr, const fortran_int_t ldvr, fortran_int_t& ilo,
  100. fortran_int_t& ihi, float* lscale, float* rscale, float& abnrm,
  101. float& bbnrm, float* rconde, float* rcondv, std::complex<float>* work,
  102. const fortran_int_t lwork, float* rwork, fortran_int_t* iwork,
  103. fortran_bool_t* bwork ) {
  104. fortran_int_t info(0);
  105. LAPACK_CGGEVX( &balanc, &jobvl, &jobvr, &sense, &n, a, &lda, b, &ldb,
  106. alpha, beta, vl, &ldvl, vr, &ldvr, &ilo, &ihi, lscale, rscale,
  107. &abnrm, &bbnrm, rconde, rcondv, work, &lwork, rwork, iwork, bwork,
  108. &info );
  109. return info;
  110. }
  111. //
  112. // Overloaded function for dispatching to
  113. // * netlib-compatible LAPACK backend (the default), and
  114. // * complex<double> value-type.
  115. //
  116. inline std::ptrdiff_t ggevx( const char balanc, const char jobvl,
  117. const char jobvr, const char sense, const fortran_int_t n,
  118. std::complex<double>* a, const fortran_int_t lda,
  119. std::complex<double>* b, const fortran_int_t ldb,
  120. std::complex<double>* alpha, std::complex<double>* beta,
  121. std::complex<double>* vl, const fortran_int_t ldvl,
  122. std::complex<double>* vr, const fortran_int_t ldvr,
  123. fortran_int_t& ilo, fortran_int_t& ihi, double* lscale,
  124. double* rscale, double& abnrm, double& bbnrm, double* rconde,
  125. double* rcondv, std::complex<double>* work, const fortran_int_t lwork,
  126. double* rwork, fortran_int_t* iwork, fortran_bool_t* bwork ) {
  127. fortran_int_t info(0);
  128. LAPACK_ZGGEVX( &balanc, &jobvl, &jobvr, &sense, &n, a, &lda, b, &ldb,
  129. alpha, beta, vl, &ldvl, vr, &ldvr, &ilo, &ihi, lscale, rscale,
  130. &abnrm, &bbnrm, rconde, rcondv, work, &lwork, rwork, iwork, bwork,
  131. &info );
  132. return info;
  133. }
  134. } // namespace detail
  135. //
  136. // Value-type based template class. Use this class if you need a type
  137. // for dispatching to ggevx.
  138. //
  139. template< typename Value, typename Enable = void >
  140. struct ggevx_impl {};
  141. //
  142. // This implementation is enabled if Value is a real type.
  143. //
  144. template< typename Value >
  145. struct ggevx_impl< Value, typename boost::enable_if< is_real< Value > >::type > {
  146. typedef Value value_type;
  147. typedef typename remove_imaginary< Value >::type real_type;
  148. //
  149. // Static member function for user-defined workspaces, that
  150. // * Deduces the required arguments for dispatching to LAPACK, and
  151. // * Asserts that most arguments make sense.
  152. //
  153. template< typename MatrixA, typename MatrixB, typename VectorALPHAR,
  154. typename VectorALPHAI, typename VectorBETA, typename MatrixVL,
  155. typename MatrixVR, typename VectorLSCALE, typename VectorRSCALE,
  156. typename VectorRCONDE, typename VectorRCONDV, typename WORK,
  157. typename IWORK, typename BWORK >
  158. static std::ptrdiff_t invoke( const char balanc, const char jobvl,
  159. const char jobvr, const char sense, MatrixA& a, MatrixB& b,
  160. VectorALPHAR& alphar, VectorALPHAI& alphai, VectorBETA& beta,
  161. MatrixVL& vl, MatrixVR& vr, fortran_int_t& ilo,
  162. fortran_int_t& ihi, VectorLSCALE& lscale,
  163. VectorRSCALE& rscale, real_type& abnrm, real_type& bbnrm,
  164. VectorRCONDE& rconde, VectorRCONDV& rcondv, detail::workspace3<
  165. WORK, IWORK, BWORK > work ) {
  166. namespace bindings = ::boost::numeric::bindings;
  167. BOOST_STATIC_ASSERT( (bindings::is_column_major< MatrixA >::value) );
  168. BOOST_STATIC_ASSERT( (bindings::is_column_major< MatrixB >::value) );
  169. BOOST_STATIC_ASSERT( (bindings::is_column_major< MatrixVL >::value) );
  170. BOOST_STATIC_ASSERT( (bindings::is_column_major< MatrixVR >::value) );
  171. BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
  172. typename bindings::value_type< MatrixA >::type >::type,
  173. typename remove_const< typename bindings::value_type<
  174. MatrixB >::type >::type >::value) );
  175. BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
  176. typename bindings::value_type< MatrixA >::type >::type,
  177. typename remove_const< typename bindings::value_type<
  178. VectorALPHAR >::type >::type >::value) );
  179. BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
  180. typename bindings::value_type< MatrixA >::type >::type,
  181. typename remove_const< typename bindings::value_type<
  182. VectorALPHAI >::type >::type >::value) );
  183. BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
  184. typename bindings::value_type< MatrixA >::type >::type,
  185. typename remove_const< typename bindings::value_type<
  186. VectorBETA >::type >::type >::value) );
  187. BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
  188. typename bindings::value_type< MatrixA >::type >::type,
  189. typename remove_const< typename bindings::value_type<
  190. MatrixVL >::type >::type >::value) );
  191. BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
  192. typename bindings::value_type< MatrixA >::type >::type,
  193. typename remove_const< typename bindings::value_type<
  194. MatrixVR >::type >::type >::value) );
  195. BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
  196. typename bindings::value_type< MatrixA >::type >::type,
  197. typename remove_const< typename bindings::value_type<
  198. VectorLSCALE >::type >::type >::value) );
  199. BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
  200. typename bindings::value_type< MatrixA >::type >::type,
  201. typename remove_const< typename bindings::value_type<
  202. VectorRSCALE >::type >::type >::value) );
  203. BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
  204. typename bindings::value_type< MatrixA >::type >::type,
  205. typename remove_const< typename bindings::value_type<
  206. VectorRCONDE >::type >::type >::value) );
  207. BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
  208. typename bindings::value_type< MatrixA >::type >::type,
  209. typename remove_const< typename bindings::value_type<
  210. VectorRCONDV >::type >::type >::value) );
  211. BOOST_STATIC_ASSERT( (bindings::is_mutable< MatrixA >::value) );
  212. BOOST_STATIC_ASSERT( (bindings::is_mutable< MatrixB >::value) );
  213. BOOST_STATIC_ASSERT( (bindings::is_mutable< VectorALPHAR >::value) );
  214. BOOST_STATIC_ASSERT( (bindings::is_mutable< VectorALPHAI >::value) );
  215. BOOST_STATIC_ASSERT( (bindings::is_mutable< VectorBETA >::value) );
  216. BOOST_STATIC_ASSERT( (bindings::is_mutable< MatrixVL >::value) );
  217. BOOST_STATIC_ASSERT( (bindings::is_mutable< MatrixVR >::value) );
  218. BOOST_STATIC_ASSERT( (bindings::is_mutable< VectorLSCALE >::value) );
  219. BOOST_STATIC_ASSERT( (bindings::is_mutable< VectorRSCALE >::value) );
  220. BOOST_STATIC_ASSERT( (bindings::is_mutable< VectorRCONDE >::value) );
  221. BOOST_STATIC_ASSERT( (bindings::is_mutable< VectorRCONDV >::value) );
  222. BOOST_ASSERT( bindings::size(alphai) >= bindings::size_column(a) );
  223. BOOST_ASSERT( bindings::size(alphar) >= bindings::size_column(a) );
  224. BOOST_ASSERT( bindings::size(work.select(fortran_int_t())) >=
  225. min_size_iwork( sense, bindings::size_column(a) ));
  226. BOOST_ASSERT( bindings::size(work.select(fortran_bool_t())) >=
  227. min_size_bwork( sense, bindings::size_column(a) ));
  228. BOOST_ASSERT( bindings::size(work.select(real_type())) >=
  229. min_size_work( balanc, jobvl, jobvr, sense,
  230. bindings::size_column(a) ));
  231. BOOST_ASSERT( bindings::size_column(a) >= 0 );
  232. BOOST_ASSERT( bindings::size_minor(a) == 1 ||
  233. bindings::stride_minor(a) == 1 );
  234. BOOST_ASSERT( bindings::size_minor(b) == 1 ||
  235. bindings::stride_minor(b) == 1 );
  236. BOOST_ASSERT( bindings::size_minor(vl) == 1 ||
  237. bindings::stride_minor(vl) == 1 );
  238. BOOST_ASSERT( bindings::size_minor(vr) == 1 ||
  239. bindings::stride_minor(vr) == 1 );
  240. BOOST_ASSERT( bindings::stride_major(a) >= std::max< std::ptrdiff_t >(1,
  241. bindings::size_column(a)) );
  242. BOOST_ASSERT( bindings::stride_major(b) >= std::max< std::ptrdiff_t >(1,
  243. bindings::size_column(a)) );
  244. BOOST_ASSERT( balanc == 'N' || balanc == 'P' || balanc == 'S' ||
  245. balanc == 'B' );
  246. BOOST_ASSERT( jobvl == 'N' || jobvl == 'V' );
  247. BOOST_ASSERT( jobvr == 'N' || jobvr == 'V' );
  248. BOOST_ASSERT( sense == 'N' || sense == 'E' || sense == 'V' ||
  249. sense == 'B' );
  250. return detail::ggevx( balanc, jobvl, jobvr, sense,
  251. bindings::size_column(a), bindings::begin_value(a),
  252. bindings::stride_major(a), bindings::begin_value(b),
  253. bindings::stride_major(b), bindings::begin_value(alphar),
  254. bindings::begin_value(alphai), bindings::begin_value(beta),
  255. bindings::begin_value(vl), bindings::stride_major(vl),
  256. bindings::begin_value(vr), bindings::stride_major(vr), ilo,
  257. ihi, bindings::begin_value(lscale),
  258. bindings::begin_value(rscale), abnrm, bbnrm,
  259. bindings::begin_value(rconde), bindings::begin_value(rcondv),
  260. bindings::begin_value(work.select(real_type())),
  261. bindings::size(work.select(real_type())),
  262. bindings::begin_value(work.select(fortran_int_t())),
  263. bindings::begin_value(work.select(fortran_bool_t())) );
  264. }
  265. //
  266. // Static member function that
  267. // * Figures out the minimal workspace requirements, and passes
  268. // the results to the user-defined workspace overload of the
  269. // invoke static member function
  270. // * Enables the unblocked algorithm (BLAS level 2)
  271. //
  272. template< typename MatrixA, typename MatrixB, typename VectorALPHAR,
  273. typename VectorALPHAI, typename VectorBETA, typename MatrixVL,
  274. typename MatrixVR, typename VectorLSCALE, typename VectorRSCALE,
  275. typename VectorRCONDE, typename VectorRCONDV >
  276. static std::ptrdiff_t invoke( const char balanc, const char jobvl,
  277. const char jobvr, const char sense, MatrixA& a, MatrixB& b,
  278. VectorALPHAR& alphar, VectorALPHAI& alphai, VectorBETA& beta,
  279. MatrixVL& vl, MatrixVR& vr, fortran_int_t& ilo,
  280. fortran_int_t& ihi, VectorLSCALE& lscale,
  281. VectorRSCALE& rscale, real_type& abnrm, real_type& bbnrm,
  282. VectorRCONDE& rconde, VectorRCONDV& rcondv, minimal_workspace ) {
  283. namespace bindings = ::boost::numeric::bindings;
  284. bindings::detail::array< real_type > tmp_work( min_size_work( balanc,
  285. jobvl, jobvr, sense, bindings::size_column(a) ) );
  286. bindings::detail::array< fortran_int_t > tmp_iwork(
  287. min_size_iwork( sense, bindings::size_column(a) ) );
  288. bindings::detail::array< fortran_bool_t > tmp_bwork( min_size_bwork(
  289. sense, bindings::size_column(a) ) );
  290. return invoke( balanc, jobvl, jobvr, sense, a, b, alphar, alphai,
  291. beta, vl, vr, ilo, ihi, lscale, rscale, abnrm, bbnrm, rconde,
  292. rcondv, workspace( tmp_work, tmp_iwork, tmp_bwork ) );
  293. }
  294. //
  295. // Static member function that
  296. // * Figures out the optimal workspace requirements, and passes
  297. // the results to the user-defined workspace overload of the
  298. // invoke static member
  299. // * Enables the blocked algorithm (BLAS level 3)
  300. //
  301. template< typename MatrixA, typename MatrixB, typename VectorALPHAR,
  302. typename VectorALPHAI, typename VectorBETA, typename MatrixVL,
  303. typename MatrixVR, typename VectorLSCALE, typename VectorRSCALE,
  304. typename VectorRCONDE, typename VectorRCONDV >
  305. static std::ptrdiff_t invoke( const char balanc, const char jobvl,
  306. const char jobvr, const char sense, MatrixA& a, MatrixB& b,
  307. VectorALPHAR& alphar, VectorALPHAI& alphai, VectorBETA& beta,
  308. MatrixVL& vl, MatrixVR& vr, fortran_int_t& ilo,
  309. fortran_int_t& ihi, VectorLSCALE& lscale,
  310. VectorRSCALE& rscale, real_type& abnrm, real_type& bbnrm,
  311. VectorRCONDE& rconde, VectorRCONDV& rcondv, optimal_workspace ) {
  312. namespace bindings = ::boost::numeric::bindings;
  313. real_type opt_size_work;
  314. bindings::detail::array< fortran_int_t > tmp_iwork(
  315. min_size_iwork( sense, bindings::size_column(a) ) );
  316. bindings::detail::array< fortran_bool_t > tmp_bwork( min_size_bwork(
  317. sense, bindings::size_column(a) ) );
  318. detail::ggevx( balanc, jobvl, jobvr, sense,
  319. bindings::size_column(a), bindings::begin_value(a),
  320. bindings::stride_major(a), bindings::begin_value(b),
  321. bindings::stride_major(b), bindings::begin_value(alphar),
  322. bindings::begin_value(alphai), bindings::begin_value(beta),
  323. bindings::begin_value(vl), bindings::stride_major(vl),
  324. bindings::begin_value(vr), bindings::stride_major(vr), ilo,
  325. ihi, bindings::begin_value(lscale),
  326. bindings::begin_value(rscale), abnrm, bbnrm,
  327. bindings::begin_value(rconde), bindings::begin_value(rcondv),
  328. &opt_size_work, -1, bindings::begin_value(tmp_iwork),
  329. bindings::begin_value(tmp_bwork) );
  330. bindings::detail::array< real_type > tmp_work(
  331. traits::detail::to_int( opt_size_work ) );
  332. return invoke( balanc, jobvl, jobvr, sense, a, b, alphar, alphai,
  333. beta, vl, vr, ilo, ihi, lscale, rscale, abnrm, bbnrm, rconde,
  334. rcondv, workspace( tmp_work, tmp_iwork, tmp_bwork ) );
  335. }
  336. //
  337. // Static member function that returns the minimum size of
  338. // workspace-array work.
  339. //
  340. static std::ptrdiff_t min_size_work( const char balanc, const char jobvl,
  341. const char jobvr, const char sense, const std::ptrdiff_t n ) {
  342. if ( balanc == 'S' || balanc == 'B' || jobvl == 'V' || jobvr == 'V' )
  343. return std::max< std::ptrdiff_t >( 1, 6*n );
  344. if ( sense == 'E' )
  345. return std::max< std::ptrdiff_t >( 1, 10*n );
  346. if ( sense == 'V' || sense == 'B' )
  347. return 2*n*n + 8*n + 16;
  348. return std::max< std::ptrdiff_t >( 1, 2*n );
  349. }
  350. //
  351. // Static member function that returns the minimum size of
  352. // workspace-array iwork.
  353. //
  354. static std::ptrdiff_t min_size_iwork( const char sense,
  355. const std::ptrdiff_t n ) {
  356. if ( sense == 'E' )
  357. return 0;
  358. else
  359. return n+6;
  360. }
  361. //
  362. // Static member function that returns the minimum size of
  363. // workspace-array bwork.
  364. //
  365. static std::ptrdiff_t min_size_bwork( const char sense,
  366. const std::ptrdiff_t n ) {
  367. if ( sense == 'N' )
  368. return 0;
  369. else
  370. return n;
  371. }
  372. };
  373. //
  374. // This implementation is enabled if Value is a complex type.
  375. //
  376. template< typename Value >
  377. struct ggevx_impl< Value, typename boost::enable_if< is_complex< Value > >::type > {
  378. typedef Value value_type;
  379. typedef typename remove_imaginary< Value >::type real_type;
  380. //
  381. // Static member function for user-defined workspaces, that
  382. // * Deduces the required arguments for dispatching to LAPACK, and
  383. // * Asserts that most arguments make sense.
  384. //
  385. template< typename MatrixA, typename MatrixB, typename VectorALPHA,
  386. typename VectorBETA, typename MatrixVL, typename MatrixVR,
  387. typename VectorLSCALE, typename VectorRSCALE,
  388. typename VectorRCONDE, typename VectorRCONDV, typename WORK,
  389. typename RWORK, typename IWORK, typename BWORK >
  390. static std::ptrdiff_t invoke( const char balanc, const char jobvl,
  391. const char jobvr, const char sense, MatrixA& a, MatrixB& b,
  392. VectorALPHA& alpha, VectorBETA& beta, MatrixVL& vl, MatrixVR& vr,
  393. fortran_int_t& ilo, fortran_int_t& ihi,
  394. VectorLSCALE& lscale, VectorRSCALE& rscale, real_type& abnrm,
  395. real_type& bbnrm, VectorRCONDE& rconde, VectorRCONDV& rcondv,
  396. detail::workspace4< WORK, RWORK, IWORK, BWORK > work ) {
  397. namespace bindings = ::boost::numeric::bindings;
  398. BOOST_STATIC_ASSERT( (bindings::is_column_major< MatrixA >::value) );
  399. BOOST_STATIC_ASSERT( (bindings::is_column_major< MatrixB >::value) );
  400. BOOST_STATIC_ASSERT( (bindings::is_column_major< MatrixVL >::value) );
  401. BOOST_STATIC_ASSERT( (bindings::is_column_major< MatrixVR >::value) );
  402. BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
  403. typename bindings::value_type< VectorLSCALE >::type >::type,
  404. typename remove_const< typename bindings::value_type<
  405. VectorRSCALE >::type >::type >::value) );
  406. BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
  407. typename bindings::value_type< VectorLSCALE >::type >::type,
  408. typename remove_const< typename bindings::value_type<
  409. VectorRCONDE >::type >::type >::value) );
  410. BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
  411. typename bindings::value_type< VectorLSCALE >::type >::type,
  412. typename remove_const< typename bindings::value_type<
  413. VectorRCONDV >::type >::type >::value) );
  414. BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
  415. typename bindings::value_type< MatrixA >::type >::type,
  416. typename remove_const< typename bindings::value_type<
  417. MatrixB >::type >::type >::value) );
  418. BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
  419. typename bindings::value_type< MatrixA >::type >::type,
  420. typename remove_const< typename bindings::value_type<
  421. VectorALPHA >::type >::type >::value) );
  422. BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
  423. typename bindings::value_type< MatrixA >::type >::type,
  424. typename remove_const< typename bindings::value_type<
  425. VectorBETA >::type >::type >::value) );
  426. BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
  427. typename bindings::value_type< MatrixA >::type >::type,
  428. typename remove_const< typename bindings::value_type<
  429. MatrixVL >::type >::type >::value) );
  430. BOOST_STATIC_ASSERT( (boost::is_same< typename remove_const<
  431. typename bindings::value_type< MatrixA >::type >::type,
  432. typename remove_const< typename bindings::value_type<
  433. MatrixVR >::type >::type >::value) );
  434. BOOST_STATIC_ASSERT( (bindings::is_mutable< MatrixA >::value) );
  435. BOOST_STATIC_ASSERT( (bindings::is_mutable< MatrixB >::value) );
  436. BOOST_STATIC_ASSERT( (bindings::is_mutable< VectorALPHA >::value) );
  437. BOOST_STATIC_ASSERT( (bindings::is_mutable< VectorBETA >::value) );
  438. BOOST_STATIC_ASSERT( (bindings::is_mutable< MatrixVL >::value) );
  439. BOOST_STATIC_ASSERT( (bindings::is_mutable< MatrixVR >::value) );
  440. BOOST_STATIC_ASSERT( (bindings::is_mutable< VectorLSCALE >::value) );
  441. BOOST_STATIC_ASSERT( (bindings::is_mutable< VectorRSCALE >::value) );
  442. BOOST_STATIC_ASSERT( (bindings::is_mutable< VectorRCONDE >::value) );
  443. BOOST_STATIC_ASSERT( (bindings::is_mutable< VectorRCONDV >::value) );
  444. BOOST_ASSERT( bindings::size(alpha) >= bindings::size_column(a) );
  445. BOOST_ASSERT( bindings::size(beta) >= bindings::size_column(a) );
  446. BOOST_ASSERT( bindings::size(work.select(fortran_int_t())) >=
  447. min_size_iwork( sense, bindings::size_column(a) ));
  448. BOOST_ASSERT( bindings::size(work.select(fortran_bool_t())) >=
  449. min_size_bwork( sense, bindings::size_column(a) ));
  450. BOOST_ASSERT( bindings::size(work.select(real_type())) >=
  451. min_size_rwork( balanc, bindings::size_column(a) ));
  452. BOOST_ASSERT( bindings::size(work.select(value_type())) >=
  453. min_size_work( sense, bindings::size_column(a) ));
  454. BOOST_ASSERT( bindings::size_column(a) >= 0 );
  455. BOOST_ASSERT( bindings::size_minor(a) == 1 ||
  456. bindings::stride_minor(a) == 1 );
  457. BOOST_ASSERT( bindings::size_minor(b) == 1 ||
  458. bindings::stride_minor(b) == 1 );
  459. BOOST_ASSERT( bindings::size_minor(vl) == 1 ||
  460. bindings::stride_minor(vl) == 1 );
  461. BOOST_ASSERT( bindings::size_minor(vr) == 1 ||
  462. bindings::stride_minor(vr) == 1 );
  463. BOOST_ASSERT( bindings::stride_major(a) >= std::max< std::ptrdiff_t >(1,
  464. bindings::size_column(a)) );
  465. BOOST_ASSERT( bindings::stride_major(b) >= std::max< std::ptrdiff_t >(1,
  466. bindings::size_column(a)) );
  467. BOOST_ASSERT( balanc == 'N' || balanc == 'P' || balanc == 'S' ||
  468. balanc == 'B' );
  469. BOOST_ASSERT( jobvl == 'N' || jobvl == 'V' );
  470. BOOST_ASSERT( jobvr == 'N' || jobvr == 'V' );
  471. BOOST_ASSERT( sense == 'N' || sense == 'E' || sense == 'V' ||
  472. sense == 'B' );
  473. return detail::ggevx( balanc, jobvl, jobvr, sense,
  474. bindings::size_column(a), bindings::begin_value(a),
  475. bindings::stride_major(a), bindings::begin_value(b),
  476. bindings::stride_major(b), bindings::begin_value(alpha),
  477. bindings::begin_value(beta), bindings::begin_value(vl),
  478. bindings::stride_major(vl), bindings::begin_value(vr),
  479. bindings::stride_major(vr), ilo, ihi,
  480. bindings::begin_value(lscale), bindings::begin_value(rscale),
  481. abnrm, bbnrm, bindings::begin_value(rconde),
  482. bindings::begin_value(rcondv),
  483. bindings::begin_value(work.select(value_type())),
  484. bindings::size(work.select(value_type())),
  485. bindings::begin_value(work.select(real_type())),
  486. bindings::begin_value(work.select(fortran_int_t())),
  487. bindings::begin_value(work.select(fortran_bool_t())) );
  488. }
  489. //
  490. // Static member function that
  491. // * Figures out the minimal workspace requirements, and passes
  492. // the results to the user-defined workspace overload of the
  493. // invoke static member function
  494. // * Enables the unblocked algorithm (BLAS level 2)
  495. //
  496. template< typename MatrixA, typename MatrixB, typename VectorALPHA,
  497. typename VectorBETA, typename MatrixVL, typename MatrixVR,
  498. typename VectorLSCALE, typename VectorRSCALE,
  499. typename VectorRCONDE, typename VectorRCONDV >
  500. static std::ptrdiff_t invoke( const char balanc, const char jobvl,
  501. const char jobvr, const char sense, MatrixA& a, MatrixB& b,
  502. VectorALPHA& alpha, VectorBETA& beta, MatrixVL& vl, MatrixVR& vr,
  503. fortran_int_t& ilo, fortran_int_t& ihi,
  504. VectorLSCALE& lscale, VectorRSCALE& rscale, real_type& abnrm,
  505. real_type& bbnrm, VectorRCONDE& rconde, VectorRCONDV& rcondv,
  506. minimal_workspace ) {
  507. namespace bindings = ::boost::numeric::bindings;
  508. bindings::detail::array< value_type > tmp_work( min_size_work( sense,
  509. bindings::size_column(a) ) );
  510. bindings::detail::array< real_type > tmp_rwork( min_size_rwork(
  511. balanc, bindings::size_column(a) ) );
  512. bindings::detail::array< fortran_int_t > tmp_iwork(
  513. min_size_iwork( sense, bindings::size_column(a) ) );
  514. bindings::detail::array< fortran_bool_t > tmp_bwork( min_size_bwork(
  515. sense, bindings::size_column(a) ) );
  516. return invoke( balanc, jobvl, jobvr, sense, a, b, alpha, beta, vl, vr,
  517. ilo, ihi, lscale, rscale, abnrm, bbnrm, rconde, rcondv,
  518. workspace( tmp_work, tmp_rwork, tmp_iwork, tmp_bwork ) );
  519. }
  520. //
  521. // Static member function that
  522. // * Figures out the optimal workspace requirements, and passes
  523. // the results to the user-defined workspace overload of the
  524. // invoke static member
  525. // * Enables the blocked algorithm (BLAS level 3)
  526. //
  527. template< typename MatrixA, typename MatrixB, typename VectorALPHA,
  528. typename VectorBETA, typename MatrixVL, typename MatrixVR,
  529. typename VectorLSCALE, typename VectorRSCALE,
  530. typename VectorRCONDE, typename VectorRCONDV >
  531. static std::ptrdiff_t invoke( const char balanc, const char jobvl,
  532. const char jobvr, const char sense, MatrixA& a, MatrixB& b,
  533. VectorALPHA& alpha, VectorBETA& beta, MatrixVL& vl, MatrixVR& vr,
  534. fortran_int_t& ilo, fortran_int_t& ihi,
  535. VectorLSCALE& lscale, VectorRSCALE& rscale, real_type& abnrm,
  536. real_type& bbnrm, VectorRCONDE& rconde, VectorRCONDV& rcondv,
  537. optimal_workspace ) {
  538. namespace bindings = ::boost::numeric::bindings;
  539. value_type opt_size_work;
  540. bindings::detail::array< real_type > tmp_rwork( min_size_rwork(
  541. balanc, bindings::size_column(a) ) );
  542. bindings::detail::array< fortran_int_t > tmp_iwork(
  543. min_size_iwork( sense, bindings::size_column(a) ) );
  544. bindings::detail::array< fortran_bool_t > tmp_bwork( min_size_bwork(
  545. sense, bindings::size_column(a) ) );
  546. detail::ggevx( balanc, jobvl, jobvr, sense,
  547. bindings::size_column(a), bindings::begin_value(a),
  548. bindings::stride_major(a), bindings::begin_value(b),
  549. bindings::stride_major(b), bindings::begin_value(alpha),
  550. bindings::begin_value(beta), bindings::begin_value(vl),
  551. bindings::stride_major(vl), bindings::begin_value(vr),
  552. bindings::stride_major(vr), ilo, ihi,
  553. bindings::begin_value(lscale), bindings::begin_value(rscale),
  554. abnrm, bbnrm, bindings::begin_value(rconde),
  555. bindings::begin_value(rcondv), &opt_size_work, -1,
  556. bindings::begin_value(tmp_rwork),
  557. bindings::begin_value(tmp_iwork),
  558. bindings::begin_value(tmp_bwork) );
  559. bindings::detail::array< value_type > tmp_work(
  560. traits::detail::to_int( opt_size_work ) );
  561. return invoke( balanc, jobvl, jobvr, sense, a, b, alpha, beta, vl, vr,
  562. ilo, ihi, lscale, rscale, abnrm, bbnrm, rconde, rcondv,
  563. workspace( tmp_work, tmp_rwork, tmp_iwork, tmp_bwork ) );
  564. }
  565. //
  566. // Static member function that returns the minimum size of
  567. // workspace-array work.
  568. //
  569. static std::ptrdiff_t min_size_work( const char sense,
  570. const std::ptrdiff_t n ) {
  571. if ( sense == 'N' )
  572. return std::max< std::ptrdiff_t >( 1, 2*n );
  573. else {
  574. if ( sense == 'E' )
  575. return std::max< std::ptrdiff_t >( 1, 4*n );
  576. else
  577. return std::max< std::ptrdiff_t >( 1, 2*n*n+2*n );
  578. }
  579. }
  580. //
  581. // Static member function that returns the minimum size of
  582. // workspace-array rwork.
  583. //
  584. static std::ptrdiff_t min_size_rwork( const char balanc,
  585. const std::ptrdiff_t n ) {
  586. if ( balanc == 'S' || balanc == 'B' )
  587. return std::max< std::ptrdiff_t >( 1, 6*n );
  588. else
  589. return std::max< std::ptrdiff_t >( 1, 2*n );
  590. }
  591. //
  592. // Static member function that returns the minimum size of
  593. // workspace-array iwork.
  594. //
  595. static std::ptrdiff_t min_size_iwork( const char sense,
  596. const std::ptrdiff_t n ) {
  597. if ( sense == 'E' )
  598. return 0;
  599. else
  600. return n+2;
  601. }
  602. //
  603. // Static member function that returns the minimum size of
  604. // workspace-array bwork.
  605. //
  606. static std::ptrdiff_t min_size_bwork( const char sense,
  607. const std::ptrdiff_t n ) {
  608. if ( sense == 'N' )
  609. return 0;
  610. else
  611. return n;
  612. }
  613. };
  614. //
  615. // Functions for direct use. These functions are overloaded for temporaries,
  616. // so that wrapped types can still be passed and used for write-access. In
  617. // addition, if applicable, they are overloaded for user-defined workspaces.
  618. // Calls to these functions are passed to the ggevx_impl classes. In the
  619. // documentation, most overloads are collapsed to avoid a large number of
  620. // prototypes which are very similar.
  621. //
  622. //
  623. // Overloaded function for ggevx. Its overload differs for
  624. // * User-defined workspace
  625. //
  626. template< typename MatrixA, typename MatrixB, typename VectorALPHAR,
  627. typename VectorALPHAI, typename VectorBETA, typename MatrixVL,
  628. typename MatrixVR, typename VectorLSCALE, typename VectorRSCALE,
  629. typename VectorRCONDE, typename VectorRCONDV, typename Workspace >
  630. inline typename boost::enable_if< detail::is_workspace< Workspace >,
  631. std::ptrdiff_t >::type
  632. ggevx( const char balanc, const char jobvl, const char jobvr,
  633. const char sense, MatrixA& a, MatrixB& b, VectorALPHAR& alphar,
  634. VectorALPHAI& alphai, VectorBETA& beta, MatrixVL& vl, MatrixVR& vr,
  635. fortran_int_t& ilo, fortran_int_t& ihi, VectorLSCALE& lscale,
  636. VectorRSCALE& rscale, typename remove_imaginary<
  637. typename bindings::value_type< MatrixA >::type >::type& abnrm,
  638. typename remove_imaginary< typename bindings::value_type<
  639. MatrixA >::type >::type& bbnrm, VectorRCONDE& rconde,
  640. VectorRCONDV& rcondv, Workspace work ) {
  641. return ggevx_impl< typename bindings::value_type<
  642. MatrixA >::type >::invoke( balanc, jobvl, jobvr, sense, a, b,
  643. alphar, alphai, beta, vl, vr, ilo, ihi, lscale, rscale, abnrm,
  644. bbnrm, rconde, rcondv, work );
  645. }
  646. //
  647. // Overloaded function for ggevx. Its overload differs for
  648. // * Default workspace-type (optimal)
  649. //
  650. template< typename MatrixA, typename MatrixB, typename VectorALPHAR,
  651. typename VectorALPHAI, typename VectorBETA, typename MatrixVL,
  652. typename MatrixVR, typename VectorLSCALE, typename VectorRSCALE,
  653. typename VectorRCONDE, typename VectorRCONDV >
  654. inline typename boost::disable_if< detail::is_workspace< VectorRCONDV >,
  655. std::ptrdiff_t >::type
  656. ggevx( const char balanc, const char jobvl, const char jobvr,
  657. const char sense, MatrixA& a, MatrixB& b, VectorALPHAR& alphar,
  658. VectorALPHAI& alphai, VectorBETA& beta, MatrixVL& vl, MatrixVR& vr,
  659. fortran_int_t& ilo, fortran_int_t& ihi, VectorLSCALE& lscale,
  660. VectorRSCALE& rscale, typename remove_imaginary<
  661. typename bindings::value_type< MatrixA >::type >::type& abnrm,
  662. typename remove_imaginary< typename bindings::value_type<
  663. MatrixA >::type >::type& bbnrm, VectorRCONDE& rconde,
  664. VectorRCONDV& rcondv ) {
  665. return ggevx_impl< typename bindings::value_type<
  666. MatrixA >::type >::invoke( balanc, jobvl, jobvr, sense, a, b,
  667. alphar, alphai, beta, vl, vr, ilo, ihi, lscale, rscale, abnrm,
  668. bbnrm, rconde, rcondv, optimal_workspace() );
  669. }
  670. //
  671. // Overloaded function for ggevx. Its overload differs for
  672. // * User-defined workspace
  673. //
  674. template< typename MatrixA, typename MatrixB, typename VectorALPHA,
  675. typename VectorBETA, typename MatrixVL, typename MatrixVR,
  676. typename VectorLSCALE, typename VectorRSCALE, typename VectorRCONDE,
  677. typename VectorRCONDV, typename Workspace >
  678. inline typename boost::enable_if< detail::is_workspace< Workspace >,
  679. std::ptrdiff_t >::type
  680. ggevx( const char balanc, const char jobvl, const char jobvr,
  681. const char sense, MatrixA& a, MatrixB& b, VectorALPHA& alpha,
  682. VectorBETA& beta, MatrixVL& vl, MatrixVR& vr, fortran_int_t& ilo,
  683. fortran_int_t& ihi, VectorLSCALE& lscale, VectorRSCALE& rscale,
  684. typename remove_imaginary< typename bindings::value_type<
  685. MatrixA >::type >::type& abnrm, typename remove_imaginary<
  686. typename bindings::value_type< MatrixA >::type >::type& bbnrm,
  687. VectorRCONDE& rconde, VectorRCONDV& rcondv, Workspace work ) {
  688. return ggevx_impl< typename bindings::value_type<
  689. MatrixA >::type >::invoke( balanc, jobvl, jobvr, sense, a, b,
  690. alpha, beta, vl, vr, ilo, ihi, lscale, rscale, abnrm, bbnrm,
  691. rconde, rcondv, work );
  692. }
  693. //
  694. // Overloaded function for ggevx. Its overload differs for
  695. // * Default workspace-type (optimal)
  696. //
  697. template< typename MatrixA, typename MatrixB, typename VectorALPHA,
  698. typename VectorBETA, typename MatrixVL, typename MatrixVR,
  699. typename VectorLSCALE, typename VectorRSCALE, typename VectorRCONDE,
  700. typename VectorRCONDV >
  701. inline typename boost::disable_if< detail::is_workspace< VectorRCONDV >,
  702. std::ptrdiff_t >::type
  703. ggevx( const char balanc, const char jobvl, const char jobvr,
  704. const char sense, MatrixA& a, MatrixB& b, VectorALPHA& alpha,
  705. VectorBETA& beta, MatrixVL& vl, MatrixVR& vr, fortran_int_t& ilo,
  706. fortran_int_t& ihi, VectorLSCALE& lscale, VectorRSCALE& rscale,
  707. typename remove_imaginary< typename bindings::value_type<
  708. MatrixA >::type >::type& abnrm, typename remove_imaginary<
  709. typename bindings::value_type< MatrixA >::type >::type& bbnrm,
  710. VectorRCONDE& rconde, VectorRCONDV& rcondv ) {
  711. return ggevx_impl< typename bindings::value_type<
  712. MatrixA >::type >::invoke( balanc, jobvl, jobvr, sense, a, b,
  713. alpha, beta, vl, vr, ilo, ihi, lscale, rscale, abnrm, bbnrm,
  714. rconde, rcondv, optimal_workspace() );
  715. }
  716. } // namespace lapack
  717. } // namespace bindings
  718. } // namespace numeric
  719. } // namespace boost
  720. #endif