[Boost-commit] svn:boost r46291 - sandbox/boost/numeric/bindings/lapack

From: karl.meerbergen_at_[hidden]
Date: 2008-06-10 07:14:53

Author: karlmeerbergen
Date: 2008-06-10 07:14:52 EDT (Tue, 10 Jun 2008)
New Revision: 46291
URL: http://svn.boost.org/trac/boost/changeset/46291

Log:
added hseqr

Added:
   sandbox/boost/numeric/bindings/lapack/hseqr.hpp   (contents, props changed)
Text files modified: 
   sandbox/boost/numeric/bindings/lapack/geqrf.hpp      |     1 +                                       
   sandbox/boost/numeric/bindings/lapack/gesvd.hpp      |     2 +-                                      
   sandbox/boost/numeric/bindings/lapack/lapack.h       |    17 +++++++++++++++++                       
   sandbox/boost/numeric/bindings/lapack/lapack_names.h |     8 ++++++++                                
   4 files changed, 27 insertions(+), 1 deletions(-)

Modified: sandbox/boost/numeric/bindings/lapack/geqrf.hpp
==============================================================================
--- sandbox/boost/numeric/bindings/lapack/geqrf.hpp	(original)
+++ sandbox/boost/numeric/bindings/lapack/geqrf.hpp	2008-06-10 07:14:52 EDT (Tue, 10 Jun 2008)
@@ -25,6 +25,7 @@
 #  include <boost/static_assert.hpp>
 #  include <boost/type_traits.hpp>
 #endif 
+#include <cassert>
 
 
 namespace boost { namespace numeric { namespace bindings { 

Modified: sandbox/boost/numeric/bindings/lapack/gesvd.hpp
==============================================================================
--- sandbox/boost/numeric/bindings/lapack/gesvd.hpp	(original)
+++ sandbox/boost/numeric/bindings/lapack/gesvd.hpp	2008-06-10 07:14:52 EDT (Tue, 10 Jun 2008)
@@ -417,7 +417,7 @@
               || (jobu == 'O' && traits::leading_dimension (u) >= 1)
               || (jobu == 'A' && traits::leading_dimension (u) >= m)
               || (jobu == 'S' && traits::leading_dimension (u) >= m));
-      assert (n == traits::matrix_size2 (vt)); 
+      assert ((jobvt == 'N' || traits::matrix_size2(vt) == n)) ;
       assert ((jobvt == 'N' && traits::leading_dimension (vt) >= 1)
               || (jobvt == 'O' && traits::leading_dimension (vt) >= 1)
               || (jobvt == 'A' && traits::leading_dimension (vt) >= n)

Added: sandbox/boost/numeric/bindings/lapack/hseqr.hpp
==============================================================================
--- (empty file)
+++ sandbox/boost/numeric/bindings/lapack/hseqr.hpp	2008-06-10 07:14:52 EDT (Tue, 10 Jun 2008)
@@ -0,0 +1,301 @@
+/*
+ * 
+ * Copyright Jeremy Conlin 2008
+ *
+ * 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 <iostream>
+#include <complex>
+#include <boost/numeric/bindings/traits/traits.hpp>
+#include <boost/numeric/bindings/traits/matrix_traits.hpp>
+#include <boost/numeric/bindings/traits/type_traits.hpp>
+#include <boost/numeric/bindings/lapack/lapack.h>
+#include <boost/numeric/bindings/traits/detail/array.hpp>
+
+#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
+#  include <boost/static_assert.hpp>
+#  include <boost/type_traits.hpp>
+#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.
+     * 
+     * There are two forms of the hseqr function:
+     *
+     * int hseqr( const char job, A& H, W& w)
+     * int hseqr( const char job, const char compz, A& H, W& w, Z& z)
+     *
+     * The first form does not compute Schur vectors and is equivelant to
+     * setting compz = 'N' in the second form.  hseqr returns a '0' if the
+     * computation is successful.
+     *
+     * job.
+     *   = 'E': compute eigenvalues only
+     *   = 'S': compute eigenvalues and the Schur form U
+     *
+     * compz. (input)
+     *   = 'N':  no Schur vectors are computed;  Equivalent to using the
+     *           first form of the hseqr function.
+     *   = '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 (output) contains the computed eigenvalues of H which is the diagonal 
+     * of U. Must be a complex object.
+     *
+     * 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 {
+        // float
+        inline
+        int hseqr_backend(const char* job, const char* compz, const int* n, 
+                const int ilo, const int ihi, float* H, const int ldH, 
+                float* wr, float* wi, float* Z, const int ldz, float* work,
+                const int* lwork){
+            int info;
+//          std::cout << "I'm inside lapack::detail::hseqr_backend for floats" 
+//              << std::endl;
+            LAPACK_SHSEQR(job, compz, n, &ilo, &ihi, H, &ldH, wr, wi, 
+                    Z, &ldz, work, lwork, &info);
+            return info;
+        }
+
+        // double
+        inline
+        int hseqr_backend(const char* job, const char* compz, const 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;
+        }
+
+        // complex<float>
+        inline
+        int hseqr_backend(const char* job, const char* compz, int* n, 
+                const int ilo, const int ihi, traits::complex_f* H, const int ldH, 
+                traits::complex_f* w, traits::complex_f* Z, int ldz, 
+                traits::complex_f* work, const int* lwork){
+            int info;
+//          std::cout << "I'm inside lapack::detail::hseqr_backend for complex<float>" 
+//              << std::endl;
+            LAPACK_CHSEQR(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;
+        }
+
+        // complex<double>
+        inline
+        int hseqr_backend(const char* job, const char* compz, int* n, 
+                const int ilo, const int ihi, traits::complex_d* H, const int ldH, 
+                traits::complex_d* w, traits::complex_d* Z, int ldz, 
+                traits::complex_d* work, const int* lwork){
+            int info;
+//          std::cout << "I'm inside lapack::detail::hseqr_backend for complex<double>" 
+//              << std::endl;
+            LAPACK_ZHSEQR(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;
+        }
+
+        template <int N>
+        struct Hseqr{};
+
+        template <>
+        struct Hseqr< 1 >{
+            template < typename A, typename W, typename V>
+            int operator() ( const char job, const char compz, A& H, W& w, V& Z ){
+//              std::cout << "Inside Hseqr<1>." << std::endl;
+
+                int n = traits::matrix_size1(H);
+                typedef typename A::value_type value_type;
+                traits::detail::array<value_type> wr(n);
+                traits::detail::array<value_type> 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<value_type> work(lwork);
+                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.storage(), &lwork);
+
+                for (int i = 0; i < n; i++)
+                    w[i] = std::complex<value_type>(wr[i], wi[i]);
+
+                return result;
+            }
+        };
+
+        template <>
+        struct Hseqr< 2 >{
+            template < typename A, typename W, typename V>
+            int operator() ( const char job, const char compz, A& H, W& w, V& Z ){
+//              std::cout << "Inside Hseqr<2>." << std::endl;
+
+                int n = traits::matrix_size1(H);
+                typedef typename A::value_type value_type;
+
+                // 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), 
+                        traits::vector_storage(w),
+                        traits::matrix_storage(Z), traits::leading_dimension(Z),
+                        &work_temp, &lwork);
+
+                if( result !=0 ) return result;
+
+                lwork = (int) std::real(work_temp);
+                traits::detail::array<value_type> work(lwork);
+                result = detail::hseqr_backend(&job, &compz, &n, 1, n,
+                        traits::matrix_storage(H),
+                        traits::leading_dimension(H), 
+                        traits::vector_storage(w),
+                        traits::matrix_storage(Z), traits::leading_dimension(Z),
+                        work.storage(), &lwork);
+
+                return result;
+            }
+        };
+        
+        template < typename A, typename W, typename V>
+        int hseqr( const char job, const char compz, A& H, W& w, V& Z ){
+//          std::cout << "I'm inside lapack::detail::hseqr." << std::endl;
+
+            assert ( job == 'E' || job == 'S' );
+            assert ( compz == 'N' || compz == 'I' || compz == 'V' );
+
+            typedef typename A::value_type value_type;
+
+            int result = detail::Hseqr< n_workspace_args<value_type>::value >()(
+                    job, compz, H, w, Z);
+
+            return result;
+        }
+    }   // namespace detail 
+
+    // Compute eigenvalues without the Schur vectors
+    template < typename A, typename W>
+    int hseqr( const char job, A& H, W& w){
+      // input checking
+#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
+      BOOST_STATIC_ASSERT((boost::is_same<
+			   typename traits::matrix_traits<A>::matrix_structure, 
+			   traits::general_t
+			   >::value)); 
+#endif 
+
+#ifndef NDEBUG
+        int const n = traits::matrix_size1(H);
+#endif
+
+        typedef typename A::value_type value_type;
+        typedef typename W::value_type complex_value_type;
+
+        assert(traits::matrix_size2(H) == n); // Square matrix
+        assert(traits::vector_size(w) == n);  
+
+        ublas::matrix<value_type, ublas::column_major> Z(1,1);
+        return detail::hseqr( job, 'N', H, w, Z );
+    }
+
+    // Compute eigenvalues and the Schur vectors
+    template < typename A, typename W, typename Z>
+    int hseqr( const char job, const char compz, A& H, W& w, Z& z){
+      // input checking
+#ifndef BOOST_NUMERIC_BINDINGS_NO_STRUCTURE_CHECK 
+      BOOST_STATIC_ASSERT((boost::is_same<
+			   typename traits::matrix_traits<A>::matrix_structure, 
+			   traits::general_t
+			   >::value)); 
+#endif 
+
+#ifndef NDEBUG
+        int const n = traits::matrix_size1(H);
+#endif
+
+        typedef typename A::value_type value_type;
+        assert(traits::matrix_size2(H) == n); // Square matrix
+        assert(traits::vector_size(w) == n);  
+        assert(traits::matrix_size2(z) == n);
+
+        return detail::hseqr( job, compz, H, w, z );
+    }
+
+  }
+}}}
+
+#endif

Modified: sandbox/boost/numeric/bindings/lapack/lapack.h
==============================================================================
--- sandbox/boost/numeric/bindings/lapack/lapack.h	(original)
+++ sandbox/boost/numeric/bindings/lapack/lapack.h	2008-06-10 07:14:52 EDT (Tue, 10 Jun 2008)
@@ -407,6 +407,23 @@
                      dcomplex_t* t, const int * ldt, dcomplex_t* q, const int* ldq,
                      int* ifst, const int * ilst, int* info );
 
+  /* Hessenberg matrices */
+
+  void LAPACK_SHSEQR( const char* JOB, const char* COMPZ, const int* N, const int* ILO, const int* IHI, float* H, 
+                      const int* LDH, float* WR, float* WI, float* Z, int const* LDZ,
+                      float* WORK, const int* LWORK, int* INFO ) ;
+
+  void LAPACK_CHSEQR( const char* JOB, const char* COMPZ, const int* N, const int* ILO, const int* IHI, fcomplex_t* H, 
+                      const int* LDH, fcomplex_t* W, fcomplex_t* Z, int const* LDZ,
+                      fcomplex_t* WORK, const int* LWORK, int* INFO ) ;
+
+  void LAPACK_DHSEQR( const char* JOB, const char* COMPZ, const int* N, const int* ILO, const int* IHI, double* H, 
+                      const int* LDH, double* WR, double* WI, double* Z, int const* LDZ,
+                      double* WORK, const int* LWORK, int* INFO ) ;
+
+  void LAPACK_ZHSEQR( const char* JOB, const char* COMPZ, const int* N, const int* ILO, const int* IHI, dcomplex_t* H, 
+                      const int* LDH, dcomplex_t* W, dcomplex_t* Z, int const* LDZ,
+                      dcomplex_t* WORK, const int* LWORK, int* INFO ) ;
 
   /* Hermitian tridiagonal matrices */
   

Modified: sandbox/boost/numeric/bindings/lapack/lapack_names.h
==============================================================================
--- sandbox/boost/numeric/bindings/lapack/lapack_names.h	(original)
+++ sandbox/boost/numeric/bindings/lapack/lapack_names.h	2008-06-10 07:14:52 EDT (Tue, 10 Jun 2008)
@@ -180,6 +180,14 @@
 
 
 /********************************************/
+/* eigenproblems for Hessenberg matrices */ 
+
+#define LAPACK_SHSEQR FORTRAN_ID( shseqr )
+#define LAPACK_DHSEQR FORTRAN_ID( dhseqr )
+#define LAPACK_CHSEQR FORTRAN_ID( chseqr )
+#define LAPACK_ZHSEQR FORTRAN_ID( zhseqr )
+
+/********************************************/
 /* eigenproblems for banded matrices */ 
 
 #define LAPACK_SSBEV FORTRAN_ID( ssbev )