[Boost-commit] svn:boost r58462 - sandbox/numeric_bindings/libs/numeric/bindings/lapack/test

Subject: [Boost-commit] svn:boost r58462 - sandbox/numeric_bindings/libs/numeric/bindings/lapack/test
From: thomas.klimpel_at_[hidden]
Date: 2009-12-18 17:41:37

Author: klimpel
Date: 2009-12-18 17:41:37 EST (Fri, 18 Dec 2009)
New Revision: 58462
URL: http://svn.boost.org/trac/boost/changeset/58462

Log:
Continue of regression tests update
Text files modified: 
   sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_hbev.cpp  |    37 +++++++++++++++++++++++++---            
   sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_hbevx.cpp |    52 ++++++++++++++++++++++++++++++++++++--- 
   sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_heevd.cpp |    40 +++++++++++++++++++++++++-----          
   sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_heevx.cpp |    47 +++++++++++++++++++++++++++++++++---    
   4 files changed, 157 insertions(+), 19 deletions(-)

Modified: sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_hbev.cpp
==============================================================================
--- sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_hbev.cpp	(original)
+++ sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_hbev.cpp	2009-12-18 17:41:37 EST (Fri, 18 Dec 2009)
@@ -8,7 +8,8 @@
 
 #include "ublas_heev.hpp"
 
-#include <boost/numeric/bindings/lapack/hbev.hpp>
+#include <boost/numeric/bindings/lapack/driver/hbev.hpp>
+#include <boost/numeric/bindings/lapack/driver/sbev.hpp>
 #include <boost/numeric/ublas/matrix_proxy.hpp>
 #include <boost/numeric/ublas/vector_proxy.hpp>
 #include <boost/numeric/bindings/traits/ublas_matrix.hpp>
@@ -17,12 +18,35 @@
 #include <boost/numeric/bindings/traits/ublas_hermitian.hpp>
 #include <boost/numeric/ublas/hermitian.hpp>
 #include <boost/numeric/ublas/io.hpp>
+#include <boost/type_traits/is_complex.hpp>
+#include <boost/mpl/if.hpp>
 
 #include <iostream>
 
 
 namespace ublas = boost::numeric::ublas;
 namespace lapack = boost::numeric::bindings::lapack;
+namespace traits = boost::numeric::bindings::traits;
+
+struct apply_real {
+  template< typename MatrixAB, typename VectorW, typename MatrixZ,
+        typename Workspace >
+  static inline integer_t hbev( const char jobz, const integer_t n,
+        const integer_t kd, MatrixAB& ab, VectorW& w, MatrixZ& z,
+        Workspace work ) {
+    return lapack::sbev( jobz, n, kd, ab, w, z, work );
+  }
+};
+
+struct apply_complex {
+  template< typename MatrixAB, typename VectorW, typename MatrixZ,
+        typename Workspace >
+  static inline integer_t hbev( const char jobz, const integer_t n,
+        const integer_t kd, MatrixAB& ab, VectorW& w, MatrixZ& z,
+        Workspace work ) {
+    return lapack::hbev( jobz, n, kd, ab, w, z, work );
+  }
+};
 
 
 template <typename U>
@@ -45,6 +69,7 @@
 
 template <typename T, typename W, typename UPLO, typename Orientation>
 int do_memory_uplo(int n, W& workspace ) {
+   typedef typename boost::mpl::if_<boost::is_complex<T>, apply_complex, apply_real>::type apply_t;
    typedef typename boost::numeric::bindings::traits::type_traits<T>::real_type real_type ;
 
    typedef ublas::banded_matrix<T, Orientation>      banded_type ;
@@ -63,11 +88,14 @@
    ublas::hermitian_adaptor<banded_type, UPLO> h( a ), h2( a2 );
 
    // Compute Schur decomposition.
-   lapack::hbev( h, e1, z, workspace ) ;
+   apply_t::hbev( 'V', traits::matrix_size2( h ), traits::matrix_upper_bandwidth( h ),
+     h, e1, z, workspace ) ;
 
    if (check_residual( h2, e1, z )) return 255 ;
 
-   lapack::hbev( h2, e2, workspace ) ;
+   matrix_type dummy_z( n, n );
+   apply_t::hbev( 'N', traits::matrix_size2( h2 ), traits::matrix_upper_bandwidth( h2 ),
+     h2, e2, dummy_z, workspace ) ;
    if (norm_2( e1 - e2 ) > n * norm_2( e1 ) * std::numeric_limits< real_type >::epsilon()) return 255 ;
 
    // Test for a matrix range
@@ -81,7 +109,8 @@
    ublas::vector_range< vector_type> e_r( e1, r );
    ublas::matrix_range< matrix_type> z_r( z, r, r );
 
-   lapack::hbev( h_r, e_r, z_r, workspace );
+   apply_t::hbev( 'V', traits::matrix_size2( h_r ), traits::matrix_upper_bandwidth( h_r ),
+     h_r, e_r, z_r, workspace );
 
    banded_range a2_r( a2, r, r );
    ublas::hermitian_adaptor< banded_range, UPLO> h2_r( a2_r );

Modified: sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_hbevx.cpp
==============================================================================
--- sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_hbevx.cpp	(original)
+++ sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_hbevx.cpp	2009-12-18 17:41:37 EST (Fri, 18 Dec 2009)
@@ -9,7 +9,8 @@
 
 #include "ublas_heev.hpp"
 
-#include <boost/numeric/bindings/lapack/hbevx.hpp>
+#include <boost/numeric/bindings/lapack/driver/hbevx.hpp>
+#include <boost/numeric/bindings/lapack/driver/sbevx.hpp>
 #include <boost/numeric/ublas/matrix_proxy.hpp>
 #include <boost/numeric/ublas/vector_proxy.hpp>
 #include <boost/numeric/bindings/traits/ublas_matrix.hpp>
@@ -18,12 +19,51 @@
 #include <boost/numeric/bindings/traits/ublas_hermitian.hpp>
 #include <boost/numeric/ublas/hermitian.hpp>
 #include <boost/numeric/ublas/io.hpp>
+#include <boost/type_traits/is_complex.hpp>
+#include <boost/mpl/if.hpp>
 
 #include <iostream>
 
 
 namespace ublas = boost::numeric::ublas;
 namespace lapack = boost::numeric::bindings::lapack;
+namespace traits = boost::numeric::bindings::traits;
+
+struct apply_real {
+  template< typename MatrixAB, typename MatrixQ, typename VectorW,
+        typename MatrixZ, typename VectorIFAIL, typename Workspace >
+  static inline integer_t hbevx( const char jobz, const char range,
+        const integer_t n, const integer_t kd, MatrixAB& ab, MatrixQ& q,
+        const typename traits::type_traits< typename traits::matrix_traits<
+        MatrixAB >::value_type >::real_type vl,
+        const typename traits::type_traits< typename traits::matrix_traits<
+        MatrixAB >::value_type >::real_type vu, const integer_t il,
+        const integer_t iu, const typename traits::type_traits<
+        typename traits::matrix_traits<
+        MatrixAB >::value_type >::real_type abstol, integer_t& m, VectorW& w,
+        MatrixZ& z, VectorIFAIL& ifail, Workspace work ) {
+    return lapack::sbevx( jobz, range, n, kd, ab, q, vl, vu,
+            il, iu, abstol, m, w, z, ifail, work );
+  }
+};
+
+struct apply_complex {
+  template< typename MatrixAB, typename MatrixQ, typename VectorW,
+        typename MatrixZ, typename VectorIFAIL, typename Workspace >
+  static inline integer_t hbevx( const char jobz, const char range,
+        const integer_t n, const integer_t kd, MatrixAB& ab, MatrixQ& q,
+        const typename traits::type_traits< typename traits::matrix_traits<
+        MatrixAB >::value_type >::real_type vl,
+        const typename traits::type_traits< typename traits::matrix_traits<
+        MatrixAB >::value_type >::real_type vu, const integer_t il,
+        const integer_t iu, const typename traits::type_traits<
+        typename traits::matrix_traits<
+        MatrixAB >::value_type >::real_type abstol, integer_t& m, VectorW& w,
+        MatrixZ& z, VectorIFAIL& ifail, Workspace work ) {
+    return lapack::hbevx( jobz, range, n, kd, ab, q, vl, vu,
+            il, iu, abstol, m, w, z, ifail, work );
+  }
+};
 
 
 template <typename U>
@@ -46,6 +86,7 @@
 
 template <typename T, typename W, typename UPLO, typename Orientation>
 int do_memory_uplo(int n, W& workspace ) {
+   typedef typename boost::mpl::if_<boost::is_complex<T>, apply_complex, apply_real>::type apply_t;
    typedef typename boost::numeric::bindings::traits::type_traits<T>::real_type real_type ;
 
    typedef ublas::banded_matrix<T, Orientation>      banded_type ;
@@ -69,11 +110,13 @@
 
 
    // Compute Schur decomposition.
-   lapack::hbevx( 'V', 'A', h, q, vl, vu, il, iu, abstol, m, e1, z, ifail, workspace ) ;
+   apply_t::hbevx( 'V', 'A', traits::matrix_size2( h ), traits::matrix_upper_bandwidth( h ),
+     h, q, vl, vu, il, iu, abstol, m, e1, z, ifail, workspace ) ;
 
    if (check_residual( h2, e1, z )) return 255 ;
 
-   lapack::hbevx( 'N', 'A', h2, q, vl, vu, il, iu, abstol, m, e2, z, ifail, workspace ) ;
+   apply_t::hbevx( 'N', 'A', traits::matrix_size2( h2 ), traits::matrix_upper_bandwidth( h2 ),
+     h2, q, vl, vu, il, iu, abstol, 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
@@ -88,7 +131,8 @@
    ublas::matrix_range< matrix_type> z_r( z, r, r );
    ublas::vector<integer_t> ifail_r(n-2);
 
-   lapack::hbevx( 'V', 'A', h_r, q, vl, vu, il, iu, abstol, m, e_r, z_r, ifail_r, workspace ) ;
+   apply_t::hbevx( 'V', 'A', traits::matrix_size2( h_r ), traits::matrix_upper_bandwidth( h_r ),
+     h_r, q, vl, vu, il, iu, abstol, m, e_r, z_r, ifail_r, workspace ) ;
 
    banded_range a2_r( a2, r, r );
    ublas::hermitian_adaptor< banded_range, UPLO> h2_r( a2_r );

Modified: sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_heevd.cpp
==============================================================================
--- sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_heevd.cpp	(original)
+++ sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_heevd.cpp	2009-12-18 17:41:37 EST (Fri, 18 Dec 2009)
@@ -9,12 +9,16 @@
 
 #include "ublas_heev.hpp"
 
-#include <boost/numeric/bindings/lapack/heevd.hpp>
+#include <boost/numeric/bindings/lapack/driver/heevd.hpp>
+#include <boost/numeric/bindings/lapack/driver/syevd.hpp>
 #include <boost/numeric/bindings/traits/ublas_matrix.hpp>
+#include <boost/numeric/bindings/traits/ublas_hermitian.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 <boost/type_traits/is_complex.hpp>
+#include <boost/mpl/if.hpp>
 
 #include <iostream>
 
@@ -22,12 +26,30 @@
 namespace ublas = boost::numeric::ublas;
 namespace lapack = boost::numeric::bindings::lapack;
 
+struct apply_real {
+  template< typename MatrixA, typename VectorW, typename Workspace >
+  static inline integer_t heevd( const char jobz, MatrixA& a, VectorW& w,
+        Workspace work ) {
+    return lapack::syevd( jobz, a, w, work );
+  }
+};
+
+struct apply_complex {
+  template< typename MatrixA, typename VectorW, typename Workspace >
+  static inline integer_t heevd( const char jobz, MatrixA& a, VectorW& w,
+        Workspace work ) {
+    return lapack::heevd( jobz, a, w, work );
+  }
+};
+
 
-template <typename T, typename W, char UPLO>
+template <typename T, typename W, typename UPLO>
 int do_memory_uplo(int n, W& workspace ) {
+   typedef typename boost::mpl::if_<boost::is_complex<T>, apply_complex, apply_real>::type apply_t;
    typedef typename boost::numeric::bindings::traits::type_traits<T>::real_type real_type ;
 
    typedef ublas::matrix<T, ublas::column_major>     matrix_type ;
+   typedef ublas::hermitian_adaptor<matrix_type, UPLO> hermitian_type ;
    typedef ublas::vector<real_type>                  vector_type ;
 
    // Set matrix
@@ -41,23 +63,27 @@
 
    // Compute Schur decomposition.
    
-   lapack::heevd( 'V', UPLO, a, e1, workspace ) ;
+   hermitian_type h_a( a );
+   apply_t::heevd( 'V', h_a, e1, workspace ) ;
 
    if (check_residual( a2, e1, a )) return 255 ;
 
-   lapack::heevd( 'N', UPLO, a2, e2, workspace ) ;
+   hermitian_type h_a2( a2 );
+   apply_t::heevd( 'N', h_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 ;
+   typedef ublas::hermitian_adaptor<matrix_range, UPLO> hermitian_range_type ;
 
    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 );
+   hermitian_range_type h_a_r( a_r );
+   apply_t::heevd( 'V', h_a_r, e_r, workspace );
 
    matrix_range a2_r( a2, r, r );
    if (check_residual( a2_r, e_r, a_r )) return 255 ;
@@ -69,9 +95,9 @@
 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 ;
+   if (do_memory_uplo<T,W,ublas::upper>(n, workspace)) return 255 ;
    std::cout << "  lower\n" ;
-   if (do_memory_uplo<T,W,'L'>(n, workspace)) return 255 ;
+   if (do_memory_uplo<T,W,ublas::lower>(n, workspace)) return 255 ;
    return 0 ;
 }
 

Modified: sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_heevx.cpp
==============================================================================
--- sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_heevx.cpp	(original)
+++ sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_heevx.cpp	2009-12-18 17:41:37 EST (Fri, 18 Dec 2009)
@@ -9,23 +9,62 @@
 
 #include "ublas_heev.hpp"
 
-#include <boost/numeric/bindings/lapack/heevx.hpp>
+#include <boost/numeric/bindings/lapack/driver/heevx.hpp>
+#include <boost/numeric/bindings/lapack/driver/syevx.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 <boost/type_traits/is_complex.hpp>
+#include <boost/mpl/if.hpp>
 
 #include <iostream>
 
 
 namespace ublas = boost::numeric::ublas;
 namespace lapack = boost::numeric::bindings::lapack;
+namespace traits = boost::numeric::bindings::traits;
+
+struct apply_real {
+  template< typename MatrixA, typename VectorW, typename MatrixZ,
+        typename VectorIFAIL, typename Workspace >
+  static inline integer_t heevx( const char jobz, const char range, MatrixA& a,
+        const typename traits::type_traits< typename traits::matrix_traits<
+        MatrixA >::value_type >::real_type vl,
+        const typename traits::type_traits< typename traits::matrix_traits<
+        MatrixA >::value_type >::real_type vu, const integer_t il,
+        const integer_t iu, const typename traits::type_traits<
+        typename traits::matrix_traits<
+        MatrixA >::value_type >::real_type abstol, integer_t& m, VectorW& w,
+        MatrixZ& z, VectorIFAIL& ifail, Workspace work ) {
+    return lapack::syevx( jobz, range, a, vl, vu, il, iu,
+            abstol, m, w, z, ifail, work );
+  }
+};
+
+struct apply_complex {
+  template< typename MatrixA, typename VectorW, typename MatrixZ,
+        typename VectorIFAIL, typename Workspace >
+  static inline integer_t heevx( const char jobz, const char range, MatrixA& a,
+        const typename traits::type_traits< typename traits::matrix_traits<
+        MatrixA >::value_type >::real_type vl,
+        const typename traits::type_traits< typename traits::matrix_traits<
+        MatrixA >::value_type >::real_type vu, const integer_t il,
+        const integer_t iu, const typename traits::type_traits<
+        typename traits::matrix_traits<
+        MatrixA >::value_type >::real_type abstol, integer_t& m, VectorW& w,
+        MatrixZ& z, VectorIFAIL& ifail, Workspace work ) {
+    return lapack::heevx( jobz, range, a, vl, vu, il, iu,
+            abstol, m, w, z, ifail, work );
+  }
+};
 
 
 template <typename T, typename W, typename UPLO>
 int do_memory_uplo(int n, W& workspace ) {
+   typedef typename boost::mpl::if_<boost::is_complex<T>, apply_complex, apply_real>::type apply_t;
    typedef typename boost::numeric::bindings::traits::type_traits<T>::real_type real_type ;
 
    typedef ublas::matrix<T, ublas::column_major>     matrix_type ;
@@ -47,13 +86,13 @@
    ublas::vector<integer_t> 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,
+   apply_t::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,
+   apply_t::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 ;
 
@@ -70,7 +109,7 @@
    ublas::vector<integer_t> 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,
+   apply_t::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 );