|
|
Boost-Commit : |
From: thomas.klimpel_at_[hidden]
Date: 2008-05-15 15:06:12
Author: klimpel
Date: 2008-05-15 15:06:12 EDT (Thu, 15 May 2008)
New Revision: 45399
URL: http://svn.boost.org/trac/boost/changeset/45399
Log:
added bindings for syevd, heevd, syevx and heevx
Added:
sandbox/libs/numeric/bindings/lapack/test/ublas_heevd.cpp (contents, props changed)
sandbox/libs/numeric/bindings/lapack/test/ublas_heevx.cpp (contents, props changed)
Text files modified:
sandbox/libs/numeric/bindings/lapack/test/Jamfile.v2 | 6 ++++++
1 files changed, 6 insertions(+), 0 deletions(-)
Modified: sandbox/libs/numeric/bindings/lapack/test/Jamfile.v2
==============================================================================
--- sandbox/libs/numeric/bindings/lapack/test/Jamfile.v2 (original)
+++ sandbox/libs/numeric/bindings/lapack/test/Jamfile.v2 2008-05-15 15:06:12 EDT (Thu, 15 May 2008)
@@ -18,6 +18,9 @@
#exe ublas_heev : ublas_heev.cpp ;
#exe ublas_syev : ublas_syev.cpp ;
+#exe ublas_heevd : ublas_heevd.cpp ;
+#exe ublas_heevx : ublas_heevx.cpp ;
+
#exe ublas_gesdd : ublas_gesdd.cc ;
#exe ublas_gesdd2 : ublas_gesdd2.cc ;
#exe ublas_gesdd3 : ublas_gesdd3.cc ;
@@ -53,6 +56,9 @@
[ run ublas_heev.cpp ]
[ run ublas_syev.cpp ]
+ [ run ublas_heevd.cpp ]
+ [ run ublas_heevx.cpp ]
+
[ run ublas_gesdd.cc ]
[ run ublas_gesdd2.cc ]
[ run ublas_gesdd3.cc ]
Added: sandbox/libs/numeric/bindings/lapack/test/ublas_heevd.cpp
==============================================================================
--- (empty file)
+++ sandbox/libs/numeric/bindings/lapack/test/ublas_heevd.cpp 2008-05-15 15:06:12 EDT (Thu, 15 May 2008)
@@ -0,0 +1,203 @@
+//
+// Copyright Toon Knapen, Karl Meerbergen
+// Copyright Thomas Klimpel 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)
+//
+
+#include "../../blas/test/random.hpp"
+
+#include <boost/numeric/bindings/lapack/heevd.hpp>
+#include <boost/numeric/bindings/traits/ublas_matrix.hpp>
+#include <boost/numeric/bindings/traits/ublas_vector.hpp>
+#include <boost/numeric/ublas/matrix_proxy.hpp>
+#include <boost/numeric/ublas/vector_proxy.hpp>
+#include <boost/numeric/ublas/io.hpp>
+
+#include <iostream>
+#include <limits>
+
+
+namespace ublas = boost::numeric::ublas;
+namespace lapack = boost::numeric::bindings::lapack;
+
+inline float conj(float v) { return v; }
+inline double conj(double v) { return v; }
+
+// Fill a matrix
+template <typename M>
+void fill(M& m) {
+ typedef typename M::size_type size_type ;
+ typedef typename M::value_type value_type ;
+
+ typedef typename boost::numeric::bindings::traits::type_traits<value_type>::real_type real_type ;
+
+ int size = m.size2() ;
+
+ for (int i=0; i<size; ++i) {
+ for (int j=0; j<i; ++j) {
+ m(j,i) = random_value<value_type>();
+ m(i,j) = conj( m(j,i) ) ;
+ }
+ m(i,i) = random_value<real_type>();
+ }
+} // randomize()
+
+
+template <typename H, typename E, typename Z>
+int check_residual(H const& h, E const& e, Z const& z) {
+ typedef typename H::value_type value_type ;
+ typedef typename E::value_type real_type ;
+
+ // Check eigen decomposition
+ int n = h.size1();
+ ublas::matrix<value_type> error( n, n ); error.clear();
+
+ // Copy band matrix in error
+ error.assign( h );
+ assert( norm_frobenius( error - herm( error ) ) == 0.0 ) ;
+
+ for (int i=0; i<n; ++i) {
+ error .minus_assign( outer_prod( column(z, i), e(i) * conj( column(z, i) ) ) ) ;
+ }
+ return (norm_frobenius( error )
+ >= n* norm_2( e ) * std::numeric_limits< real_type >::epsilon() ) ;
+} // check_residual()
+
+
+template <typename T, typename W, char UPLO>
+int do_memory_uplo(int n, W& workspace ) {
+ typedef typename boost::numeric::bindings::traits::type_traits<T>::real_type real_type ;
+
+ typedef ublas::matrix<T, ublas::column_major> matrix_type ;
+ typedef ublas::vector<real_type> vector_type ;
+
+ // Set matrix
+ matrix_type a( n, n ); a.clear();
+ vector_type e1( n );
+ vector_type e2( n );
+
+ fill( a );
+ matrix_type a2( a );
+ matrix_type z( a );
+
+ // Compute Schur decomposition.
+
+ lapack::heevd( 'V', UPLO, a, e1, workspace ) ;
+
+ if (check_residual( a2, e1, a )) return 255 ;
+
+ lapack::heevd( 'N', UPLO, a2, e2, workspace ) ;
+ if (norm_2( e1 - e2 ) > n * norm_2( e1 ) * std::numeric_limits< real_type >::epsilon()) return 255 ;
+
+ // Test for a matrix range
+ fill( a ); a2.assign( a );
+
+ typedef ublas::matrix_range< matrix_type > matrix_range ;
+
+ ublas::range r(1,n-1) ;
+ matrix_range a_r( a, r, r );
+ ublas::vector_range< vector_type> e_r( e1, r );
+
+ lapack::heevd( 'V', UPLO, a_r, e_r, workspace );
+
+ matrix_range a2_r( a2, r, r );
+ if (check_residual( a2_r, e_r, a_r )) return 255 ;
+
+ return 0 ;
+} // do_memory_uplo()
+
+
+template <typename T, typename W>
+int do_memory_type(int n, W workspace) {
+ std::cout << " upper\n" ;
+ if (do_memory_uplo<T,W,'U'>(n, workspace)) return 255 ;
+ std::cout << " lower\n" ;
+ if (do_memory_uplo<T,W,'L'>(n, workspace)) return 255 ;
+ return 0 ;
+}
+
+
+template <typename T>
+struct Workspace {
+ typedef ublas::vector< T > array_type ;
+ typedef ublas::vector< int > int_array_type ;
+ typedef lapack::detail::workspace1< array_type > first_type;
+ typedef lapack::detail::workspace1< int_array_type > second_type;
+ typedef std::pair<first_type, second_type> type ;
+
+ Workspace(size_t n)
+ : work_( 1+6*n+2*n*n )
+ , iwork_( 3+5*n )
+ {}
+
+ type operator() () {
+ return type( first_type( work_ ), second_type( iwork_ ) );
+ }
+
+ array_type work_ ;
+ int_array_type iwork_ ;
+};
+
+template <typename T>
+struct Workspace< std::complex<T> > {
+ typedef ublas::vector< std::complex<T> > complex_array_type ;
+ typedef ublas::vector< T > real_array_type ;
+ typedef ublas::vector< int > int_array_type ;
+ typedef lapack::detail::workspace2< complex_array_type,real_array_type > first_type;
+ typedef lapack::detail::workspace1< int_array_type > second_type;
+ typedef std::pair<first_type, second_type> type ;
+
+ Workspace(size_t n)
+ : work_( 2*n+n*n )
+ , rwork_( 1+5*n+2*n*n )
+ , iwork_( 3+5*n )
+ {}
+
+ type operator() () {
+ return type( first_type( work_, rwork_ ), second_type( iwork_ ) );
+ }
+
+ complex_array_type work_ ;
+ real_array_type rwork_ ;
+ int_array_type iwork_ ;
+};
+
+
+template <typename T>
+int do_value_type() {
+ const int n = 8 ;
+
+ std::cout << " optimal workspace\n";
+ if (do_memory_type<T,lapack::optimal_workspace>( n, lapack::optimal_workspace() ) ) return 255 ;
+
+ std::cout << " minimal workspace\n";
+ if (do_memory_type<T,lapack::minimal_workspace>( n, lapack::minimal_workspace() ) ) return 255 ;
+
+ std::cout << " workspace array\n";
+ Workspace<T> work( n );
+ do_memory_type<T,typename Workspace<T>::type >( n, work() );
+ return 0;
+} // do_value_type()
+
+
+int main() {
+ // Run tests for different value_types
+ std::cout << "float\n" ;
+ if (do_value_type< float >()) return 255;
+
+ std::cout << "double\n" ;
+ if (do_value_type< double >()) return 255;
+
+ std::cout << "complex<float>\n" ;
+ if (do_value_type< std::complex<float> >()) return 255;
+
+ std::cout << "complex<double>\n" ;
+ if (do_value_type< std::complex<double> >()) return 255;
+
+ std::cout << "Regression test succeeded\n" ;
+ return 0;
+}
+
Added: sandbox/libs/numeric/bindings/lapack/test/ublas_heevx.cpp
==============================================================================
--- (empty file)
+++ sandbox/libs/numeric/bindings/lapack/test/ublas_heevx.cpp 2008-05-15 15:06:12 EDT (Thu, 15 May 2008)
@@ -0,0 +1,217 @@
+//
+// Copyright Toon Knapen, Karl Meerbergen
+// Copyright Thomas Klimpel 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)
+//
+
+#include "../../blas/test/random.hpp"
+
+#include <boost/numeric/bindings/lapack/heevx.hpp>
+#include <boost/numeric/bindings/traits/ublas_matrix.hpp>
+#include <boost/numeric/bindings/traits/ublas_vector.hpp>
+#include <boost/numeric/bindings/traits/ublas_hermitian.hpp>
+#include <boost/numeric/ublas/matrix_proxy.hpp>
+#include <boost/numeric/ublas/vector_proxy.hpp>
+#include <boost/numeric/ublas/io.hpp>
+
+#include <iostream>
+#include <limits>
+
+
+namespace ublas = boost::numeric::ublas;
+namespace lapack = boost::numeric::bindings::lapack;
+
+inline float conj(float v) { return v; }
+inline double conj(double v) { return v; }
+
+// Fill a matrix
+template <typename M>
+void fill(M& m) {
+ typedef typename M::size_type size_type ;
+ typedef typename M::value_type value_type ;
+
+ typedef typename boost::numeric::bindings::traits::type_traits<value_type>::real_type real_type ;
+
+ int size = m.size2() ;
+
+ for (int i=0; i<size; ++i) {
+ for (int j=0; j<i; ++j) {
+ m(j,i) = random_value<value_type>();
+ m(i,j) = conj( m(j,i) ) ;
+ }
+ m(i,i) = random_value<real_type>();
+ }
+} // randomize()
+
+
+template <typename H, typename E, typename Z>
+int check_residual(H const& h, E const& e, Z const& z) {
+ typedef typename H::value_type value_type ;
+ typedef typename E::value_type real_type ;
+
+ // Check eigen decomposition
+ int n = h.size1();
+ ublas::matrix<value_type> error( n, n ); error.clear();
+
+ // Copy band matrix in error
+ error.assign( h );
+ assert( norm_frobenius( error - herm( error ) ) == 0.0 ) ;
+
+ for (int i=0; i<n; ++i) {
+ error .minus_assign( outer_prod( column(z, i), e(i) * conj( column(z, i) ) ) ) ;
+ }
+ return (norm_frobenius( error )
+ >= n* norm_2( e ) * std::numeric_limits< real_type >::epsilon() ) ;
+} // check_residual()
+
+
+template <typename T, typename W, typename UPLO>
+int do_memory_uplo(int n, W& workspace ) {
+ typedef typename boost::numeric::bindings::traits::type_traits<T>::real_type real_type ;
+
+ typedef ublas::matrix<T, ublas::column_major> matrix_type ;
+ typedef ublas::vector<real_type> vector_type ;
+
+ typedef ublas::hermitian_adaptor<matrix_type, UPLO> hermitian_type;
+
+ // Set matrix
+ matrix_type a( n, n ); a.clear();
+ vector_type e1( n );
+ vector_type e2( n );
+
+ fill( a );
+ matrix_type a2( a );
+ matrix_type z( a );
+
+ // Compute Schur decomposition.
+ int m;
+ ublas::vector<int> ifail(n);
+
+ hermitian_type h_a( a );
+ lapack::heevx( 'V', 'A', h_a, real_type(0.0), real_type(1.0), 2, n-1, real_type(1e-28), m,
+ e1, z, ifail, workspace ) ;
+
+ if (check_residual( a2, e1, z )) return 255 ;
+
+ hermitian_type h_a2( a2 );
+ lapack::heevx( 'N', 'A', h_a2, real_type(0.0), real_type(1.0), 2, n-1, real_type(1e-28), m,
+ e2, z, ifail, workspace ) ;
+ if (norm_2( e1 - e2 ) > n * norm_2( e1 ) * std::numeric_limits< real_type >::epsilon()) return 255 ;
+
+ // Test for a matrix range
+ fill( a ); a2.assign( a );
+
+ typedef ublas::matrix_range< matrix_type > matrix_range ;
+ typedef ublas::hermitian_adaptor<matrix_range, UPLO> hermitian_range_type;
+
+ ublas::range r(1,n-1) ;
+ matrix_range a_r( a, r, r );
+ matrix_range z_r( z, r, r );
+ ublas::vector_range< vector_type> e_r( e1, r );
+ ublas::vector<int> ifail_r(n-2);
+
+ hermitian_range_type h_a_r( a_r );
+ lapack::heevx( 'V', 'A', h_a_r, real_type(0.0), real_type(1.0), 2, n-1, real_type(1e-28), m,
+ e_r, z_r, ifail_r, workspace );
+
+ matrix_range a2_r( a2, r, r );
+ if (check_residual( a2_r, e_r, z_r )) return 255 ;
+
+ return 0 ;
+} // do_memory_uplo()
+
+
+template <typename T, typename W>
+int do_memory_type(int n, W workspace) {
+ std::cout << " upper\n" ;
+ if (do_memory_uplo<T,W,ublas::upper>(n, workspace)) return 255 ;
+ std::cout << " lower\n" ;
+ if (do_memory_uplo<T,W,ublas::lower>(n, workspace)) return 255 ;
+ return 0 ;
+}
+
+
+template <typename T>
+struct Workspace {
+ typedef ublas::vector< T > array_type ;
+ typedef ublas::vector< int > int_array_type ;
+ typedef lapack::detail::workspace1< array_type > first_type;
+ typedef lapack::detail::workspace1< int_array_type > second_type;
+ typedef std::pair<first_type, second_type> type ;
+
+ Workspace(size_t n)
+ : work_( 8*n )
+ , iwork_( 5*n )
+ {}
+
+ type operator() () {
+ return type( first_type( work_ ), second_type( iwork_ ) );
+ }
+
+ array_type work_ ;
+ int_array_type iwork_ ;
+};
+
+template <typename T>
+struct Workspace< std::complex<T> > {
+ typedef ublas::vector< std::complex<T> > complex_array_type ;
+ typedef ublas::vector< T > real_array_type ;
+ typedef ublas::vector< int > int_array_type ;
+ typedef lapack::detail::workspace2< complex_array_type,real_array_type > first_type;
+ typedef lapack::detail::workspace1< int_array_type > second_type;
+ typedef std::pair<first_type, second_type> type ;
+
+ Workspace(size_t n)
+ : work_( 2*n )
+ , rwork_( 7*n )
+ , iwork_( 5*n )
+ {}
+
+ type operator() () {
+ return type( first_type( work_, rwork_ ), second_type( iwork_ ) );
+ }
+
+ complex_array_type work_ ;
+ real_array_type rwork_ ;
+ int_array_type iwork_ ;
+};
+
+
+template <typename T>
+int do_value_type() {
+ const int n = 8 ;
+
+ std::cout << " optimal workspace\n";
+ if (do_memory_type<T,lapack::optimal_workspace>( n, lapack::optimal_workspace() ) ) return 255 ;
+
+ std::cout << " minimal workspace\n";
+ if (do_memory_type<T,lapack::minimal_workspace>( n, lapack::minimal_workspace() ) ) return 255 ;
+
+ std::cout << " workspace array\n";
+ Workspace<T> work( n );
+ do_memory_type<T,typename Workspace<T>::type >( n, work() );
+ return 0;
+} // do_value_type()
+
+
+int main() {
+ // Run tests for different value_types
+ std::cout << "float\n" ;
+ if (do_value_type< float >()) return 255;
+
+ std::cout << "double\n" ;
+ if (do_value_type< double >()) return 255;
+
+ std::cout << "complex<float>\n" ;
+ if (do_value_type< std::complex<float> >()) return 255;
+
+ std::cout << "complex<double>\n" ;
+ if (do_value_type< std::complex<double> >()) return 255;
+
+ std::cout << "Regression test succeeded\n" ;
+ return 0;
+}
+
Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk