PageRenderTime 57ms CodeModel.GetById 24ms RepoModel.GetById 0ms app.codeStats 0ms

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

http://github.com/quantmind/jflib
C++ Header | 559 lines | 503 code | 29 blank | 27 comment | 67 complexity | c2daf76b5d43741aa27f610889c06238 MD5 | raw file
Possible License(s): BSD-3-Clause 
    

  1. //
    2. 

  2. // Copyright (c) 2003--2009
    3. 

  3. // Toon Knapen, Karl Meerbergen, Kresimir Fresl,
    4. 

  4. // Thomas Klimpel and Rutger ter Borg
    5. 

  5. //
    6. 

  6. // Distributed under the Boost Software License, Version 1.0.
    7. 

  7. // (See accompanying file LICENSE_1_0.txt or copy at
    8. 

  8. // http://www.boost.org/LICENSE_1_0.txt)
    9. 

  9. //
    10. 

  10. // THIS FILE IS AUTOMATICALLY GENERATED
    11. 

  11. // PLEASE DO NOT EDIT!
    12. 

  12. //
    13. 

  13. 

  14. #ifndef BOOST_NUMERIC_BINDINGS_LAPACK_DRIVER_GGEVX_HPP
    15. 

  15. #define BOOST_NUMERIC_BINDINGS_LAPACK_DRIVER_GGEVX_HPP
    16. 

  16. 

  17. #include <boost/assert.hpp>
    18. 

  18. #include <boost/mpl/bool.hpp>
    19. 

  19. #include <boost/numeric/bindings/lapack/detail/lapack.h>
    20. 

  20. #include <boost/numeric/bindings/lapack/workspace.hpp>
    21. 

  21. #include <boost/numeric/bindings/traits/detail/array.hpp>
    22. 

  22. #include <boost/numeric/bindings/traits/detail/utils.hpp>
    23. 

  23. #include <boost/numeric/bindings/traits/is_complex.hpp>
    24. 

  24. #include <boost/numeric/bindings/traits/is_real.hpp>
    25. 

  25. #include <boost/numeric/bindings/traits/traits.hpp>
    26. 

  26. #include <boost/numeric/bindings/traits/type_traits.hpp>
    27. 

  27. #include <boost/static_assert.hpp>
    28. 

  28. #include <boost/type_traits/is_same.hpp>
    29. 

  29. #include <boost/utility/enable_if.hpp>
    30. 

  30. 

  31. namespace boost {
    32. 

  32. namespace numeric {
    33. 

  33. namespace bindings {
    34. 

  34. namespace lapack {
    35. 

  35. 

  36. //$DESCRIPTION
    37. 

  37. 

  38. // overloaded functions to call lapack
    39. 

  39. namespace detail {
    40. 

  40.     inline void ggevx( char const balanc, char const jobvl, char const jobvr,
    41. 

  41.             char const sense, integer_t const n, float* a,
    42. 

  42.             integer_t const lda, float* b, integer_t const ldb, float* alphar,
    43. 

  43.             float* alphai, float* beta, float* vl, integer_t const ldvl,
    44. 

  44.             float* vr, integer_t const ldvr, integer_t& ilo, integer_t& ihi,
    45. 

  45.             float* lscale, float* rscale, float& abnrm, float& bbnrm,
    46. 

  46.             float* rconde, float* rcondv, float* work, integer_t const lwork,
    47. 

  47.             integer_t* iwork, logical_t* bwork, integer_t& info ) {
    48. 

  48.         LAPACK_SGGEVX( &balanc, &jobvl, &jobvr, &sense, &n, a, &lda, b, &ldb,
    49. 

  49.                 alphar, alphai, beta, vl, &ldvl, vr, &ldvr, &ilo, &ihi,
    50. 

  50.                 lscale, rscale, &abnrm, &bbnrm, rconde, rcondv, work, &lwork,
    51. 

  51.                 iwork, bwork, &info );
    52. 

  52.     }
    53. 

  53.     inline void ggevx( char const balanc, char const jobvl, char const jobvr,
    54. 

  54.             char const sense, integer_t const n, double* a,
    55. 

  55.             integer_t const lda, double* b, integer_t const ldb,
    56. 

  56.             double* alphar, double* alphai, double* beta, double* vl,
    57. 

  57.             integer_t const ldvl, double* vr, integer_t const ldvr,
    58. 

  58.             integer_t& ilo, integer_t& ihi, double* lscale, double* rscale,
    59. 

  59.             double& abnrm, double& bbnrm, double* rconde, double* rcondv,
    60. 

  60.             double* work, integer_t const lwork, integer_t* iwork,
    61. 

  61.             logical_t* bwork, integer_t& info ) {
    62. 

  62.         LAPACK_DGGEVX( &balanc, &jobvl, &jobvr, &sense, &n, a, &lda, b, &ldb,
    63. 

  63.                 alphar, alphai, beta, vl, &ldvl, vr, &ldvr, &ilo, &ihi,
    64. 

  64.                 lscale, rscale, &abnrm, &bbnrm, rconde, rcondv, work, &lwork,
    65. 

  65.                 iwork, bwork, &info );
    66. 

  66.     }
    67. 

  67.     inline void ggevx( char const balanc, char const jobvl, char const jobvr,
    68. 

  68.             char const sense, integer_t const n, traits::complex_f* a,
    69. 

  69.             integer_t const lda, traits::complex_f* b, integer_t const ldb,
    70. 

  70.             traits::complex_f* alpha, traits::complex_f* beta,
    71. 

  71.             traits::complex_f* vl, integer_t const ldvl,
    72. 

  72.             traits::complex_f* vr, integer_t const ldvr, integer_t& ilo,
    73. 

  73.             integer_t& ihi, float* lscale, float* rscale, float& abnrm,
    74. 

  74.             float& bbnrm, float* rconde, float* rcondv,
    75. 

  75.             traits::complex_f* work, integer_t const lwork, float* rwork,
    76. 

  76.             integer_t* iwork, logical_t* bwork, integer_t& info ) {
    77. 

  77.         LAPACK_CGGEVX( &balanc, &jobvl, &jobvr, &sense, &n,
    78. 

  78.                 traits::complex_ptr(a), &lda, traits::complex_ptr(b), &ldb,
    79. 

  79.                 traits::complex_ptr(alpha), traits::complex_ptr(beta),
    80. 

  80.                 traits::complex_ptr(vl), &ldvl, traits::complex_ptr(vr),
    81. 

  81.                 &ldvr, &ilo, &ihi, lscale, rscale, &abnrm, &bbnrm, rconde,
    82. 

  82.                 rcondv, traits::complex_ptr(work), &lwork, rwork, iwork,
    83. 

  83.                 bwork, &info );
    84. 

  84.     }
    85. 

  85.     inline void ggevx( char const balanc, char const jobvl, char const jobvr,
    86. 

  86.             char const sense, integer_t const n, traits::complex_d* a,
    87. 

  87.             integer_t const lda, traits::complex_d* b, integer_t const ldb,
    88. 

  88.             traits::complex_d* alpha, traits::complex_d* beta,
    89. 

  89.             traits::complex_d* vl, integer_t const ldvl,
    90. 

  90.             traits::complex_d* vr, integer_t const ldvr, integer_t& ilo,
    91. 

  91.             integer_t& ihi, double* lscale, double* rscale, double& abnrm,
    92. 

  92.             double& bbnrm, double* rconde, double* rcondv,
    93. 

  93.             traits::complex_d* work, integer_t const lwork, double* rwork,
    94. 

  94.             integer_t* iwork, logical_t* bwork, integer_t& info ) {
    95. 

  95.         LAPACK_ZGGEVX( &balanc, &jobvl, &jobvr, &sense, &n,
    96. 

  96.                 traits::complex_ptr(a), &lda, traits::complex_ptr(b), &ldb,
    97. 

  97.                 traits::complex_ptr(alpha), traits::complex_ptr(beta),
    98. 

  98.                 traits::complex_ptr(vl), &ldvl, traits::complex_ptr(vr),
    99. 

  99.                 &ldvr, &ilo, &ihi, lscale, rscale, &abnrm, &bbnrm, rconde,
    100. 

  100.                 rcondv, traits::complex_ptr(work), &lwork, rwork, iwork,
    101. 

  101.                 bwork, &info );
    102. 

  102.     }
    103. 

  103. }
    104. 

  104. 

  105. // value-type based template
    106. 

  106. template< typename ValueType, typename Enable = void >
    107. 

  107. struct ggevx_impl{};
    108. 

  108. 

  109. // real specialization
    110. 

  110. template< typename ValueType >
    111. 

  111. struct ggevx_impl< ValueType, typename boost::enable_if< traits::is_real<ValueType> >::type > {
    112. 

  112. 

  113.     typedef ValueType value_type;
    114. 

  114.     typedef typename traits::type_traits<ValueType>::real_type real_type;
    115. 

  115. 

  116.     // user-defined workspace specialization
    117. 

  117.     template< typename MatrixA, typename MatrixB, typename VectorALPHAR,
    118. 

  118.             typename VectorALPHAI, typename VectorBETA, typename MatrixVL,
    119. 

  119.             typename MatrixVR, typename VectorLSCALE, typename VectorRSCALE,
    120. 

  120.             typename VectorRCONDE, typename VectorRCONDV, typename WORK,
    121. 

  121.             typename IWORK, typename BWORK >
    122. 

  122.     static void invoke( char const balanc, char const jobvl, char const jobvr,
    123. 

  123.             char const sense, MatrixA& a, MatrixB& b, VectorALPHAR& alphar,
    124. 

  124.             VectorALPHAI& alphai, VectorBETA& beta, MatrixVL& vl,
    125. 

  125.             MatrixVR& vr, integer_t& ilo, integer_t& ihi,
    126. 

  126.             VectorLSCALE& lscale, VectorRSCALE& rscale, real_type& abnrm,
    127. 

  127.             real_type& bbnrm, VectorRCONDE& rconde, VectorRCONDV& rcondv,
    128. 

  128.             integer_t& info, detail::workspace3< WORK, IWORK, BWORK > work ) {
    129. 

  129.         BOOST_STATIC_ASSERT( (boost::is_same< typename traits::matrix_traits<
    130. 

  130.                 MatrixA >::value_type, typename traits::matrix_traits<
    131. 

  131.                 MatrixB >::value_type >::value) );
    132. 

  132.         BOOST_STATIC_ASSERT( (boost::is_same< typename traits::matrix_traits<
    133. 

  133.                 MatrixA >::value_type, typename traits::vector_traits<
    134. 

  134.                 VectorALPHAR >::value_type >::value) );
    135. 

  135.         BOOST_STATIC_ASSERT( (boost::is_same< typename traits::matrix_traits<
    136. 

  136.                 MatrixA >::value_type, typename traits::vector_traits<
    137. 

  137.                 VectorALPHAI >::value_type >::value) );
    138. 

  138.         BOOST_STATIC_ASSERT( (boost::is_same< typename traits::matrix_traits<
    139. 

  139.                 MatrixA >::value_type, typename traits::vector_traits<
    140. 

  140.                 VectorBETA >::value_type >::value) );
    141. 

  141.         BOOST_STATIC_ASSERT( (boost::is_same< typename traits::matrix_traits<
    142. 

  142.                 MatrixA >::value_type, typename traits::matrix_traits<
    143. 

  143.                 MatrixVL >::value_type >::value) );
    144. 

  144.         BOOST_STATIC_ASSERT( (boost::is_same< typename traits::matrix_traits<
    145. 

  145.                 MatrixA >::value_type, typename traits::matrix_traits<
    146. 

  146.                 MatrixVR >::value_type >::value) );
    147. 

  147.         BOOST_STATIC_ASSERT( (boost::is_same< typename traits::matrix_traits<
    148. 

  148.                 MatrixA >::value_type, typename traits::vector_traits<
    149. 

  149.                 VectorLSCALE >::value_type >::value) );
    150. 

  150.         BOOST_STATIC_ASSERT( (boost::is_same< typename traits::matrix_traits<
    151. 

  151.                 MatrixA >::value_type, typename traits::vector_traits<
    152. 

  152.                 VectorRSCALE >::value_type >::value) );
    153. 

  153.         BOOST_STATIC_ASSERT( (boost::is_same< typename traits::matrix_traits<
    154. 

  154.                 MatrixA >::value_type, typename traits::vector_traits<
    155. 

  155.                 VectorRCONDE >::value_type >::value) );
    156. 

  156.         BOOST_STATIC_ASSERT( (boost::is_same< typename traits::matrix_traits<
    157. 

  157.                 MatrixA >::value_type, typename traits::vector_traits<
    158. 

  158.                 VectorRCONDV >::value_type >::value) );
    159. 

  159.         BOOST_ASSERT( balanc == 'N' || balanc == 'P' || balanc == 'S' ||
    160. 

  160.                 balanc == 'B' );
    161. 

  161.         BOOST_ASSERT( jobvl == 'N' || jobvl == 'V' );
    162. 

  162.         BOOST_ASSERT( jobvr == 'N' || jobvr == 'V' );
    163. 

  163.         BOOST_ASSERT( sense == 'N' || sense == 'E' || sense == 'V' ||
    164. 

  164.                 sense == 'B' );
    165. 

  165.         BOOST_ASSERT( traits::matrix_num_columns(a) >= 0 );
    166. 

  166.         BOOST_ASSERT( traits::leading_dimension(a) >= std::max(1,
    167. 

  167.                 traits::matrix_num_columns(a)) );
    168. 

  168.         BOOST_ASSERT( traits::leading_dimension(b) >= std::max(1,
    169. 

  169.                 traits::matrix_num_columns(a)) );
    170. 

  170.         BOOST_ASSERT( traits::vector_size(alphar) >=
    171. 

  171.                 traits::matrix_num_columns(a) );
    172. 

  172.         BOOST_ASSERT( traits::vector_size(alphai) >=
    173. 

  173.                 traits::matrix_num_columns(a) );
    174. 

  174.         BOOST_ASSERT( traits::vector_size(work.select(real_type())) >=
    175. 

  175.                 min_size_work( balanc, jobvl, jobvr, sense,
    176. 

  176.                 traits::matrix_num_columns(a) ));
    177. 

  177.         BOOST_ASSERT( traits::vector_size(work.select(integer_t())) >=
    178. 

  178.                 min_size_iwork( sense, traits::matrix_num_columns(a) ));
    179. 

  179.         BOOST_ASSERT( traits::vector_size(work.select(bool())) >=
    180. 

  180.                 min_size_bwork( sense, traits::matrix_num_columns(a) ));
    181. 

  181.         detail::ggevx( balanc, jobvl, jobvr, sense,
    182. 

  182.                 traits::matrix_num_columns(a), traits::matrix_storage(a),
    183. 

  183.                 traits::leading_dimension(a), traits::matrix_storage(b),
    184. 

  184.                 traits::leading_dimension(b), traits::vector_storage(alphar),
    185. 

  185.                 traits::vector_storage(alphai), traits::vector_storage(beta),
    186. 

  186.                 traits::matrix_storage(vl), traits::leading_dimension(vl),
    187. 

  187.                 traits::matrix_storage(vr), traits::leading_dimension(vr),
    188. 

  188.                 ilo, ihi, traits::vector_storage(lscale),
    189. 

  189.                 traits::vector_storage(rscale), abnrm, bbnrm,
    190. 

  190.                 traits::vector_storage(rconde),
    191. 

  191.                 traits::vector_storage(rcondv),
    192. 

  192.                 traits::vector_storage(work.select(real_type())),
    193. 

  193.                 traits::vector_size(work.select(real_type())),
    194. 

  194.                 traits::vector_storage(work.select(integer_t())),
    195. 

  195.                 traits::vector_storage(work.select(bool())), info );
    196. 

  196.     }
    197. 

  197. 

  198.     // minimal workspace specialization
    199. 

  199.     template< typename MatrixA, typename MatrixB, typename VectorALPHAR,
    200. 

  200.             typename VectorALPHAI, typename VectorBETA, typename MatrixVL,
    201. 

  201.             typename MatrixVR, typename VectorLSCALE, typename VectorRSCALE,
    202. 

  202.             typename VectorRCONDE, typename VectorRCONDV >
    203. 

  203.     static void invoke( char const balanc, char const jobvl, char const jobvr,
    204. 

  204.             char const sense, MatrixA& a, MatrixB& b, VectorALPHAR& alphar,
    205. 

  205.             VectorALPHAI& alphai, VectorBETA& beta, MatrixVL& vl,
    206. 

  206.             MatrixVR& vr, integer_t& ilo, integer_t& ihi,
    207. 

  207.             VectorLSCALE& lscale, VectorRSCALE& rscale, real_type& abnrm,
    208. 

  208.             real_type& bbnrm, VectorRCONDE& rconde, VectorRCONDV& rcondv,
    209. 

  209.             integer_t& info, minimal_workspace work ) {
    210. 

  210.         traits::detail::array< real_type > tmp_work( min_size_work( balanc,
    211. 

  211.                 jobvl, jobvr, sense, traits::matrix_num_columns(a) ) );
    212. 

  212.         traits::detail::array< integer_t > tmp_iwork( min_size_iwork( sense,
    213. 

  213.                 traits::matrix_num_columns(a) ) );
    214. 

  214.         traits::detail::array< bool > tmp_bwork( min_size_bwork( sense,
    215. 

  215.                 traits::matrix_num_columns(a) ) );
    216. 

  216.         invoke( balanc, jobvl, jobvr, sense, a, b, alphar, alphai, beta, vl,
    217. 

  217.                 vr, ilo, ihi, lscale, rscale, abnrm, bbnrm, rconde, rcondv,
    218. 

  218.                 info, workspace( tmp_work, tmp_iwork, tmp_bwork ) );
    219. 

  219.     }
    220. 

  220. 

  221.     // optimal workspace specialization
    222. 

  222.     template< typename MatrixA, typename MatrixB, typename VectorALPHAR,
    223. 

  223.             typename VectorALPHAI, typename VectorBETA, typename MatrixVL,
    224. 

  224.             typename MatrixVR, typename VectorLSCALE, typename VectorRSCALE,
    225. 

  225.             typename VectorRCONDE, typename VectorRCONDV >
    226. 

  226.     static void invoke( char const balanc, char const jobvl, char const jobvr,
    227. 

  227.             char const sense, MatrixA& a, MatrixB& b, VectorALPHAR& alphar,
    228. 

  228.             VectorALPHAI& alphai, VectorBETA& beta, MatrixVL& vl,
    229. 

  229.             MatrixVR& vr, integer_t& ilo, integer_t& ihi,
    230. 

  230.             VectorLSCALE& lscale, VectorRSCALE& rscale, real_type& abnrm,
    231. 

  231.             real_type& bbnrm, VectorRCONDE& rconde, VectorRCONDV& rcondv,
    232. 

  232.             integer_t& info, optimal_workspace work ) {
    233. 

  233.         real_type opt_size_work;
    234. 

  234.         traits::detail::array< integer_t > tmp_iwork( min_size_iwork( sense,
    235. 

  235.                 traits::matrix_num_columns(a) ) );
    236. 

  236.         traits::detail::array< bool > tmp_bwork( min_size_bwork( sense,
    237. 

  237.                 traits::matrix_num_columns(a) ) );
    238. 

  238.         detail::ggevx( balanc, jobvl, jobvr, sense,
    239. 

  239.                 traits::matrix_num_columns(a), traits::matrix_storage(a),
    240. 

  240.                 traits::leading_dimension(a), traits::matrix_storage(b),
    241. 

  241.                 traits::leading_dimension(b), traits::vector_storage(alphar),
    242. 

  242.                 traits::vector_storage(alphai), traits::vector_storage(beta),
    243. 

  243.                 traits::matrix_storage(vl), traits::leading_dimension(vl),
    244. 

  244.                 traits::matrix_storage(vr), traits::leading_dimension(vr),
    245. 

  245.                 ilo, ihi, traits::vector_storage(lscale),
    246. 

  246.                 traits::vector_storage(rscale), abnrm, bbnrm,
    247. 

  247.                 traits::vector_storage(rconde),
    248. 

  248.                 traits::vector_storage(rcondv), &opt_size_work, -1,
    249. 

  249.                 traits::vector_storage(tmp_iwork),
    250. 

  250.                 traits::vector_storage(tmp_bwork), info );
    251. 

  251.         traits::detail::array< real_type > tmp_work(
    252. 

  252.                 traits::detail::to_int( opt_size_work ) );
    253. 

  253.         invoke( balanc, jobvl, jobvr, sense, a, b, alphar, alphai, beta, vl,
    254. 

  254.                 vr, ilo, ihi, lscale, rscale, abnrm, bbnrm, rconde, rcondv,
    255. 

  255.                 info, workspace( tmp_work, tmp_iwork, tmp_bwork ) );
    256. 

  256.     }
    257. 

  257. 

  258.     static integer_t min_size_work( char const balanc, char const jobvl,
    259. 

  259.             char const jobvr, char const sense, integer_t const n ) {
    260. 

  260.         if ( balanc == 'S' || balanc == 'B' || jobvl == 'V' || jobvr == 'V' )
    261. 

  261.             return std::max( 1, 6*n );
    262. 

  262.         if ( sense == 'E' )
    263. 

  263.             return std::max( 1, 10*n );
    264. 

  264.         if ( sense == 'V' || sense == 'B' )
    265. 

  265.             return 2*n*n + 8*n + 16;
    266. 

  266.         return std::max( 1, 2*n );
    267. 

  267.     }
    268. 

  268. 

  269.     static integer_t min_size_iwork( char const sense, integer_t const n ) {
    270. 

  270.         if ( sense == 'E' )
    271. 

  271.           return 0;
    272. 

  272.         else
    273. 

  273.           return n+6;
    274. 

  274.     }
    275. 

  275. 

  276.     static integer_t min_size_bwork( char const sense, integer_t const n ) {
    277. 

  277.         if ( sense == 'N' )
    278. 

  278.           return 0;
    279. 

  279.         else
    280. 

  280.           return n;
    281. 

  281.     }
    282. 

  282. };
    283. 

  283. 

  284. // complex specialization
    285. 

  285. template< typename ValueType >
    286. 

  286. struct ggevx_impl< ValueType, typename boost::enable_if< traits::is_complex<ValueType> >::type > {
    287. 

  287. 

  288.     typedef ValueType value_type;
    289. 

  289.     typedef typename traits::type_traits<ValueType>::real_type real_type;
    290. 

  290. 

  291.     // user-defined workspace specialization
    292. 

  292.     template< typename MatrixA, typename MatrixB, typename VectorALPHA,
    293. 

  293.             typename VectorBETA, typename MatrixVL, typename MatrixVR,
    294. 

  294.             typename VectorLSCALE, typename VectorRSCALE,
    295. 

  295.             typename VectorRCONDE, typename VectorRCONDV, typename WORK,
    296. 

  296.             typename RWORK, typename IWORK, typename BWORK >
    297. 

  297.     static void invoke( char const balanc, char const jobvl, char const jobvr,
    298. 

  298.             char const sense, MatrixA& a, MatrixB& b, VectorALPHA& alpha,
    299. 

  299.             VectorBETA& beta, MatrixVL& vl, MatrixVR& vr, integer_t& ilo,
    300. 

  300.             integer_t& ihi, VectorLSCALE& lscale, VectorRSCALE& rscale,
    301. 

  301.             real_type& abnrm, real_type& bbnrm, VectorRCONDE& rconde,
    302. 

  302.             VectorRCONDV& rcondv, integer_t& info, detail::workspace4< WORK,
    303. 

  303.             RWORK, IWORK, BWORK > work ) {
    304. 

  304.         BOOST_STATIC_ASSERT( (boost::is_same< typename traits::vector_traits<
    305. 

  305.                 VectorLSCALE >::value_type, typename traits::vector_traits<
    306. 

  306.                 VectorRSCALE >::value_type >::value) );
    307. 

  307.         BOOST_STATIC_ASSERT( (boost::is_same< typename traits::vector_traits<
    308. 

  308.                 VectorLSCALE >::value_type, typename traits::vector_traits<
    309. 

  309.                 VectorRCONDE >::value_type >::value) );
    310. 

  310.         BOOST_STATIC_ASSERT( (boost::is_same< typename traits::vector_traits<
    311. 

  311.                 VectorLSCALE >::value_type, typename traits::vector_traits<
    312. 

  312.                 VectorRCONDV >::value_type >::value) );
    313. 

  313.         BOOST_STATIC_ASSERT( (boost::is_same< typename traits::matrix_traits<
    314. 

  314.                 MatrixA >::value_type, typename traits::matrix_traits<
    315. 

  315.                 MatrixB >::value_type >::value) );
    316. 

  316.         BOOST_STATIC_ASSERT( (boost::is_same< typename traits::matrix_traits<
    317. 

  317.                 MatrixA >::value_type, typename traits::vector_traits<
    318. 

  318.                 VectorALPHA >::value_type >::value) );
    319. 

  319.         BOOST_STATIC_ASSERT( (boost::is_same< typename traits::matrix_traits<
    320. 

  320.                 MatrixA >::value_type, typename traits::vector_traits<
    321. 

  321.                 VectorBETA >::value_type >::value) );
    322. 

  322.         BOOST_STATIC_ASSERT( (boost::is_same< typename traits::matrix_traits<
    323. 

  323.                 MatrixA >::value_type, typename traits::matrix_traits<
    324. 

  324.                 MatrixVL >::value_type >::value) );
    325. 

  325.         BOOST_STATIC_ASSERT( (boost::is_same< typename traits::matrix_traits<
    326. 

  326.                 MatrixA >::value_type, typename traits::matrix_traits<
    327. 

  327.                 MatrixVR >::value_type >::value) );
    328. 

  328.         BOOST_ASSERT( balanc == 'N' || balanc == 'P' || balanc == 'S' ||
    329. 

  329.                 balanc == 'B' );
    330. 

  330.         BOOST_ASSERT( jobvl == 'N' || jobvl == 'V' );
    331. 

  331.         BOOST_ASSERT( jobvr == 'N' || jobvr == 'V' );
    332. 

  332.         BOOST_ASSERT( sense == 'N' || sense == 'E' || sense == 'V' ||
    333. 

  333.                 sense == 'B' );
    334. 

  334.         BOOST_ASSERT( traits::matrix_num_columns(a) >= 0 );
    335. 

  335.         BOOST_ASSERT( traits::leading_dimension(a) >= std::max(1,
    336. 

  336.                 traits::matrix_num_columns(a)) );
    337. 

  337.         BOOST_ASSERT( traits::leading_dimension(b) >= std::max(1,
    338. 

  338.                 traits::matrix_num_columns(a)) );
    339. 

  339.         BOOST_ASSERT( traits::vector_size(alpha) >=
    340. 

  340.                 traits::matrix_num_columns(a) );
    341. 

  341.         BOOST_ASSERT( traits::vector_size(beta) >=
    342. 

  342.                 traits::matrix_num_columns(a) );
    343. 

  343.         BOOST_ASSERT( traits::vector_size(work.select(value_type())) >=
    344. 

  344.                 min_size_work( sense, traits::matrix_num_columns(a) ));
    345. 

  345.         BOOST_ASSERT( traits::vector_size(work.select(real_type())) >=
    346. 

  346.                 min_size_rwork( balanc, traits::matrix_num_columns(a) ));
    347. 

  347.         BOOST_ASSERT( traits::vector_size(work.select(integer_t())) >=
    348. 

  348.                 min_size_iwork( sense, traits::matrix_num_columns(a) ));
    349. 

  349.         BOOST_ASSERT( traits::vector_size(work.select(bool())) >=
    350. 

  350.                 min_size_bwork( sense, traits::matrix_num_columns(a) ));
    351. 

  351.         detail::ggevx( balanc, jobvl, jobvr, sense,
    352. 

  352.                 traits::matrix_num_columns(a), traits::matrix_storage(a),
    353. 

  353.                 traits::leading_dimension(a), traits::matrix_storage(b),
    354. 

  354.                 traits::leading_dimension(b), traits::vector_storage(alpha),
    355. 

  355.                 traits::vector_storage(beta), traits::matrix_storage(vl),
    356. 

  356.                 traits::leading_dimension(vl), traits::matrix_storage(vr),
    357. 

  357.                 traits::leading_dimension(vr), ilo, ihi,
    358. 

  358.                 traits::vector_storage(lscale),
    359. 

  359.                 traits::vector_storage(rscale), abnrm, bbnrm,
    360. 

  360.                 traits::vector_storage(rconde),
    361. 

  361.                 traits::vector_storage(rcondv),
    362. 

  362.                 traits::vector_storage(work.select(value_type())),
    363. 

  363.                 traits::vector_size(work.select(value_type())),
    364. 

  364.                 traits::vector_storage(work.select(real_type())),
    365. 

  365.                 traits::vector_storage(work.select(integer_t())),
    366. 

  366.                 traits::vector_storage(work.select(bool())), info );
    367. 

  367.     }
    368. 

  368. 

  369.     // minimal workspace specialization
    370. 

  370.     template< typename MatrixA, typename MatrixB, typename VectorALPHA,
    371. 

  371.             typename VectorBETA, typename MatrixVL, typename MatrixVR,
    372. 

  372.             typename VectorLSCALE, typename VectorRSCALE,
    373. 

  373.             typename VectorRCONDE, typename VectorRCONDV >
    374. 

  374.     static void invoke( char const balanc, char const jobvl, char const jobvr,
    375. 

  375.             char const sense, MatrixA& a, MatrixB& b, VectorALPHA& alpha,
    376. 

  376.             VectorBETA& beta, MatrixVL& vl, MatrixVR& vr, integer_t& ilo,
    377. 

  377.             integer_t& ihi, VectorLSCALE& lscale, VectorRSCALE& rscale,
    378. 

  378.             real_type& abnrm, real_type& bbnrm, VectorRCONDE& rconde,
    379. 

  379.             VectorRCONDV& rcondv, integer_t& info, minimal_workspace work ) {
    380. 

  380.         traits::detail::array< value_type > tmp_work( min_size_work( sense,
    381. 

  381.                 traits::matrix_num_columns(a) ) );
    382. 

  382.         traits::detail::array< real_type > tmp_rwork( min_size_rwork( balanc,
    383. 

  383.                 traits::matrix_num_columns(a) ) );
    384. 

  384.         traits::detail::array< integer_t > tmp_iwork( min_size_iwork( sense,
    385. 

  385.                 traits::matrix_num_columns(a) ) );
    386. 

  386.         traits::detail::array< bool > tmp_bwork( min_size_bwork( sense,
    387. 

  387.                 traits::matrix_num_columns(a) ) );
    388. 

  388.         invoke( balanc, jobvl, jobvr, sense, a, b, alpha, beta, vl, vr, ilo,
    389. 

  389.                 ihi, lscale, rscale, abnrm, bbnrm, rconde, rcondv, info,
    390. 

  390.                 workspace( tmp_work, tmp_rwork, tmp_iwork, tmp_bwork ) );
    391. 

  391.     }
    392. 

  392. 

  393.     // optimal workspace specialization
    394. 

  394.     template< typename MatrixA, typename MatrixB, typename VectorALPHA,
    395. 

  395.             typename VectorBETA, typename MatrixVL, typename MatrixVR,
    396. 

  396.             typename VectorLSCALE, typename VectorRSCALE,
    397. 

  397.             typename VectorRCONDE, typename VectorRCONDV >
    398. 

  398.     static void invoke( char const balanc, char const jobvl, char const jobvr,
    399. 

  399.             char const sense, MatrixA& a, MatrixB& b, VectorALPHA& alpha,
    400. 

  400.             VectorBETA& beta, MatrixVL& vl, MatrixVR& vr, integer_t& ilo,
    401. 

  401.             integer_t& ihi, VectorLSCALE& lscale, VectorRSCALE& rscale,
    402. 

  402.             real_type& abnrm, real_type& bbnrm, VectorRCONDE& rconde,
    403. 

  403.             VectorRCONDV& rcondv, integer_t& info, optimal_workspace work ) {
    404. 

  404.         value_type opt_size_work;
    405. 

  405.         traits::detail::array< real_type > tmp_rwork( min_size_rwork( balanc,
    406. 

  406.                 traits::matrix_num_columns(a) ) );
    407. 

  407.         traits::detail::array< integer_t > tmp_iwork( min_size_iwork( sense,
    408. 

  408.                 traits::matrix_num_columns(a) ) );
    409. 

  409.         traits::detail::array< bool > tmp_bwork( min_size_bwork( sense,
    410. 

  410.                 traits::matrix_num_columns(a) ) );
    411. 

  411.         detail::ggevx( balanc, jobvl, jobvr, sense,
    412. 

  412.                 traits::matrix_num_columns(a), traits::matrix_storage(a),
    413. 

  413.                 traits::leading_dimension(a), traits::matrix_storage(b),
    414. 

  414.                 traits::leading_dimension(b), traits::vector_storage(alpha),
    415. 

  415.                 traits::vector_storage(beta), traits::matrix_storage(vl),
    416. 

  416.                 traits::leading_dimension(vl), traits::matrix_storage(vr),
    417. 

  417.                 traits::leading_dimension(vr), ilo, ihi,
    418. 

  418.                 traits::vector_storage(lscale),
    419. 

  419.                 traits::vector_storage(rscale), abnrm, bbnrm,
    420. 

  420.                 traits::vector_storage(rconde),
    421. 

  421.                 traits::vector_storage(rcondv), &opt_size_work, -1,
    422. 

  422.                 traits::vector_storage(tmp_rwork),
    423. 

  423.                 traits::vector_storage(tmp_iwork),
    424. 

  424.                 traits::vector_storage(tmp_bwork), info );
    425. 

  425.         traits::detail::array< value_type > tmp_work(
    426. 

  426.                 traits::detail::to_int( opt_size_work ) );
    427. 

  427.         invoke( balanc, jobvl, jobvr, sense, a, b, alpha, beta, vl, vr, ilo,
    428. 

  428.                 ihi, lscale, rscale, abnrm, bbnrm, rconde, rcondv, info,
    429. 

  429.                 workspace( tmp_work, tmp_rwork, tmp_iwork, tmp_bwork ) );
    430. 

  430.     }
    431. 

  431. 

  432.     static integer_t min_size_work( char const sense, integer_t const n ) {
    433. 

  433.         if ( sense == 'N' )
    434. 

  434.             return std::max( 1, 2*n );
    435. 

  435.         else {
    436. 

  436.             if ( sense == 'E' )
    437. 

  437.                 return std::max( 1, 4*n );
    438. 

  438.             else
    439. 

  439.                 return std::max( 1, 2*n*n+2*n );
    440. 

  440.         }
    441. 

  441.     }
    442. 

  442. 

  443.     static integer_t min_size_rwork( char const balanc, integer_t const n ) {
    444. 

  444.         if ( balanc == 'S' || balanc == 'B' )
    445. 

  445.             return std::max( 1, 6*n );
    446. 

  446.         else
    447. 

  447.             return std::max( 1, 2*n );
    448. 

  448.     }
    449. 

  449. 

  450.     static integer_t min_size_iwork( char const sense, integer_t const n ) {
    451. 

  451.         if ( sense == 'E' )
    452. 

  452.           return 0;
    453. 

  453.         else
    454. 

  454.           return n+2;
    455. 

  455.     }
    456. 

  456. 

  457.     static integer_t min_size_bwork( char const sense, integer_t const n ) {
    458. 

  458.         if ( sense == 'N' )
    459. 

  459.           return 0;
    460. 

  460.         else
    461. 

  461.           return n;
    462. 

  462.     }
    463. 

  463. };
    464. 

  464. 

  465. 

  466. // template function to call ggevx
    467. 

  467. template< typename MatrixA, typename MatrixB, typename VectorALPHAR,
    468. 

  468.         typename VectorALPHAI, typename VectorBETA, typename MatrixVL,
    469. 

  469.         typename MatrixVR, typename VectorLSCALE, typename VectorRSCALE,
    470. 

  470.         typename VectorRCONDE, typename VectorRCONDV, typename Workspace >
    471. 

  471. inline integer_t ggevx( char const balanc, char const jobvl,
    472. 

  472.         char const jobvr, char const sense, MatrixA& a, MatrixB& b,
    473. 

  473.         VectorALPHAR& alphar, VectorALPHAI& alphai, VectorBETA& beta,
    474. 

  474.         MatrixVL& vl, MatrixVR& vr, integer_t& ilo, integer_t& ihi,
    475. 

  475.         VectorLSCALE& lscale, VectorRSCALE& rscale,
    476. 

  476.         typename traits::type_traits< typename traits::matrix_traits<
    477. 

  477.         MatrixA >::value_type >::real_type& abnrm,
    478. 

  478.         typename traits::type_traits< typename traits::matrix_traits<
    479. 

  479.         MatrixA >::value_type >::real_type& bbnrm, VectorRCONDE& rconde,
    480. 

  480.         VectorRCONDV& rcondv, Workspace work ) {
    481. 

  481.     typedef typename traits::matrix_traits< MatrixA >::value_type value_type;
    482. 

  482.     integer_t info(0);
    483. 

  483.     ggevx_impl< value_type >::invoke( balanc, jobvl, jobvr, sense, a, b,
    484. 

  484.             alphar, alphai, beta, vl, vr, ilo, ihi, lscale, rscale, abnrm,
    485. 

  485.             bbnrm, rconde, rcondv, info, work );
    486. 

  486.     return info;
    487. 

  487. }
    488. 

  488. 

  489. // template function to call ggevx, default workspace type
    490. 

  490. template< typename MatrixA, typename MatrixB, typename VectorALPHAR,
    491. 

  491.         typename VectorALPHAI, typename VectorBETA, typename MatrixVL,
    492. 

  492.         typename MatrixVR, typename VectorLSCALE, typename VectorRSCALE,
    493. 

  493.         typename VectorRCONDE, typename VectorRCONDV >
    494. 

  494. inline integer_t ggevx( char const balanc, char const jobvl,
    495. 

  495.         char const jobvr, char const sense, MatrixA& a, MatrixB& b,
    496. 

  496.         VectorALPHAR& alphar, VectorALPHAI& alphai, VectorBETA& beta,
    497. 

  497.         MatrixVL& vl, MatrixVR& vr, integer_t& ilo, integer_t& ihi,
    498. 

  498.         VectorLSCALE& lscale, VectorRSCALE& rscale,
    499. 

  499.         typename traits::type_traits< typename traits::matrix_traits<
    500. 

  500.         MatrixA >::value_type >::real_type& abnrm,
    501. 

  501.         typename traits::type_traits< typename traits::matrix_traits<
    502. 

  502.         MatrixA >::value_type >::real_type& bbnrm, VectorRCONDE& rconde,
    503. 

  503.         VectorRCONDV& rcondv ) {
    504. 

  504.     typedef typename traits::matrix_traits< MatrixA >::value_type value_type;
    505. 

  505.     integer_t info(0);
    506. 

  506.     ggevx_impl< value_type >::invoke( balanc, jobvl, jobvr, sense, a, b,
    507. 

  507.             alphar, alphai, beta, vl, vr, ilo, ihi, lscale, rscale, abnrm,
    508. 

  508.             bbnrm, rconde, rcondv, info, optimal_workspace() );
    509. 

  509.     return info;
    510. 

  510. }
    511. 

  511. // template function to call ggevx
    512. 

  512. template< typename MatrixA, typename MatrixB, typename VectorALPHA,
    513. 

  513.         typename VectorBETA, typename MatrixVL, typename MatrixVR,
    514. 

  514.         typename VectorLSCALE, typename VectorRSCALE, typename VectorRCONDE,
    515. 

  515.         typename VectorRCONDV, typename Workspace >
    516. 

  516. inline integer_t ggevx( char const balanc, char const jobvl,
    517. 

  517.         char const jobvr, char const sense, MatrixA& a, MatrixB& b,
    518. 

  518.         VectorALPHA& alpha, VectorBETA& beta, MatrixVL& vl, MatrixVR& vr,
    519. 

  519.         integer_t& ilo, integer_t& ihi, VectorLSCALE& lscale,
    520. 

  520.         VectorRSCALE& rscale, typename traits::type_traits<
    521. 

  521.         typename traits::matrix_traits<
    522. 

  522.         MatrixA >::value_type >::real_type& abnrm,
    523. 

  523.         typename traits::type_traits< typename traits::matrix_traits<
    524. 

  524.         MatrixA >::value_type >::real_type& bbnrm, VectorRCONDE& rconde,
    525. 

  525.         VectorRCONDV& rcondv, Workspace work ) {
    526. 

  526.     typedef typename traits::matrix_traits< MatrixA >::value_type value_type;
    527. 

  527.     integer_t info(0);
    528. 

  528.     ggevx_impl< value_type >::invoke( balanc, jobvl, jobvr, sense, a, b,
    529. 

  529.             alpha, beta, vl, vr, ilo, ihi, lscale, rscale, abnrm, bbnrm,
    530. 

  530.             rconde, rcondv, info, work );
    531. 

  531.     return info;
    532. 

  532. }
    533. 

  533. 

  534. // template function to call ggevx, default workspace type
    535. 

  535. template< typename MatrixA, typename MatrixB, typename VectorALPHA,
    536. 

  536.         typename VectorBETA, typename MatrixVL, typename MatrixVR,
    537. 

  537.         typename VectorLSCALE, typename VectorRSCALE, typename VectorRCONDE,
    538. 

  538.         typename VectorRCONDV >
    539. 

  539. inline integer_t ggevx( char const balanc, char const jobvl,
    540. 

  540.         char const jobvr, char const sense, MatrixA& a, MatrixB& b,
    541. 

  541.         VectorALPHA& alpha, VectorBETA& beta, MatrixVL& vl, MatrixVR& vr,
    542. 

  542.         integer_t& ilo, integer_t& ihi, VectorLSCALE& lscale,
    543. 

  543.         VectorRSCALE& rscale, typename traits::type_traits<
    544. 

  544.         typename traits::matrix_traits<
    545. 

  545.         MatrixA >::value_type >::real_type& abnrm,
    546. 

  546.         typename traits::type_traits< typename traits::matrix_traits<
    547. 

  547.         MatrixA >::value_type >::real_type& bbnrm, VectorRCONDE& rconde,
    548. 

  548.         VectorRCONDV& rcondv ) {
    549. 

  549.     typedef typename traits::matrix_traits< MatrixA >::value_type value_type;
    550. 

  550.     integer_t info(0);
    551. 

  551.     ggevx_impl< value_type >::invoke( balanc, jobvl, jobvr, sense, a, b,
    552. 

  552.             alpha, beta, vl, vr, ilo, ihi, lscale, rscale, abnrm, bbnrm,
    553. 

  553.             rconde, rcondv, info, optimal_workspace() );
    554. 

  554.     return info;
    555. 

  555. }
    556. 

  556. 

  557. }}}} // namespace boost::numeric::bindings::lapack
    558. 

  558. 

  559. #endif
    560. 

  560.