/* * * * Distributed under the Boost Software License, Version 1.0. * (See accompanying file LICENSE_1_0.txt or copy at * http://www.boost.org/LICENSE_1_0.txt) * */ #ifndef BOOST_NUMERIC_BINDINGS_LAPACK_HSEQR_HPP #define BOOST_NUMERIC_BINDINGS_LAPACK_HSEQR_HPP #include #include #include #include #include #include #ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK # include # include #endif namespace boost { namespace numeric { namespace bindings { namespace lapack { /////////////////////////////////////////////////////////////////// // // Compute eigenvalues of an Hessenberg matrix, H. // /////////////////////////////////////////////////////////////////// /* * hseqr() computes the eigenvalues of a Hessenberg matrix H * and, optionally, the matrices T and Z from the Schur decomposition * H = Z U Z**T, where U is an upper quasi-triangular matrix (the * Schur form), and Z is the orthogonal matrix of Schur vectors. * * Optionally Z may be postmultiplied into an input orthogonal * matrix Q so that this routine can give the Schur factorization * of a matrix A which has been reduced to the Hessenberg form H * by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*U*(QZ)**T. * * int hseqr ( char job, char compz, A& H, W& w, V* Z, optimal_workspace); * * job. (input) * = 'E': compute eigenvalues only * = 'S': compute eigenvalues and the Schur form U * * compz. (input) * = 'N': no Schur vectors are computed; * = 'I': Z is initialized to the unit matrix and the matrix Z * of Schur vectors of H is returned; * = 'V': Z must contain an orthogonal matrix Q on entry, and * the product Q*Z is returned. * * H is the Hessenberg matrix whose eigenpairs you're interested * in. (input/output) On exit, if computation is successful and * job = 'S', then H contains the * upper quasi-triangular matrix U from the Schur decomposition * (the Schur form); 2-by-2 diagonal blocks (corresponding to * complex conjugate pairs of eigenvalues) are returned in * standard form, with H(i,i) = H(i+1,i+1) and * H(i+1,i)*H(i,i+1) < 0. If computation is successful and * job = 'E', the contents of H are unspecified on exit. * * w contains the computed eigenvalues of H which is the diagonal * of U. (output) * * Z. (input/output) * If compz = 'N', Z is not referenced. * If compz = 'I', on entry Z need not be set and on exit, * if computation is successful, Z contains the orthogonal matrix Z of the Schur * vectors of H. If compz = 'V', on entry Z must contain an * N-by-N matrix Q, which is assumed to be equal to the unit * matrix . On exit, if computation is successful, Z contains Q*Z. * */ namespace detail { inline int hseqr_backend(const char* job, const char* compz, const unsigned int* n, const int ilo, const int ihi, double* H, const int ldH, double* wr, double* wi, double* Z, const int ldz, double* work, const int* lwork){ int info; std::cout << "I'm inside lapack::detail::hseqr_backend for doubles" << std::endl; LAPACK_DHSEQR(job, compz, n, &ilo, &ihi, H, &ldH, wr, wi, Z, &ldz, work, lwork, &info); return info; } // inline // int hseqr_backend(const char* job, const char* compz, const unsigned int* n, // const int ilo, const int ihi, traits::complex_d* H, const int ldH, // traits::complex_d* w, traits::complex_d* Z, const int ldz, // traits::complex_d* work, const int* lwork){ // int info; // std::cout << "I'm inside lapack::detail::hseqr_backend for complexs" // << std::endl; // LAPACK_DHSEQR(job, compz, n, &ilo, &ihi, // traits::complex_ptr(H), &ldH, // traits::complex_ptr(w), // traits::complex_ptr(Z), &ldz, // traits::complex_ptr(work), lwork, &info); // return info; // } } // namespace detail template < typename A, typename W, typename V> int hseqr( char job, char compz, A& H, W& w, V& Z){ // std::cout << "I'm inside lapack::detail::hseqr." << std::endl; unsigned int const n = traits::matrix_size1(H); typedef typename A::value_type value_type; traits::detail::array wr(n); traits::detail::array wi(n); // workspace query int lwork = -1; value_type work_temp; int result = detail::hseqr_backend(&job, &compz, &n, 1, n, traits::matrix_storage(H), traits::leading_dimension(H), wr.storage(), wi.storage(), traits::matrix_storage(Z), traits::leading_dimension(Z), &work_temp, &lwork); if( result !=0 ) return result; lwork = (int) work_temp; traits::detail::array work(lwork); int results = detail::hseqr_backend(&job, &compz, &n, 1, n, traits::matrix_storage(H), traits::leading_dimension(H), wr.storage(), wi.storage(), traits::matrix_storage(Z), traits::leading_dimension(Z), &work_temp, &lwork); for (int i = 0; i < n; i++) w[i] = std::complex(wr[i], wi[i]); return result; } } }}} #endif