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

Subject: [Boost-commit] svn:boost r58474 - in sandbox/numeric_bindings: boost/numeric/bindings/lapack/driver libs/numeric/bindings/lapack/test
From: thomas.klimpel_at_[hidden]
Date: 2009-12-21 08:47:04

Author: klimpel
Date: 2009-12-21 08:47:04 EST (Mon, 21 Dec 2009)
New Revision: 58474
URL: http://svn.boost.org/trac/boost/changeset/58474

Log:
Continue of regression tests update
Text files modified: 
   sandbox/numeric_bindings/boost/numeric/bindings/lapack/driver/gees.hpp    |     4 +-                                      
   sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_gees.cpp |    58 +++++++++++++++++++++++++++++++++------ 
   sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_geev.cpp |    11 ++++--                                  
   3 files changed, 58 insertions(+), 15 deletions(-)

Modified: sandbox/numeric_bindings/boost/numeric/bindings/lapack/driver/gees.hpp
==============================================================================
--- sandbox/numeric_bindings/boost/numeric/bindings/lapack/driver/gees.hpp	(original)
+++ sandbox/numeric_bindings/boost/numeric/bindings/lapack/driver/gees.hpp	2009-12-21 08:47:04 EST (Mon, 21 Dec 2009)
@@ -268,7 +268,7 @@
 // template function to call gees
 template< typename MatrixA, typename VectorWR, typename VectorWI,
         typename MatrixVS, typename Workspace >
-inline integer_t gees( const char jobvs, const char sort,
+inline integer_t gees_2( const char jobvs, const char sort,
         logical_t* select, MatrixA& a, integer_t& sdim, VectorWR& wr,
         VectorWI& wi, MatrixVS& vs, Workspace work ) {
     typedef typename traits::matrix_traits< MatrixA >::value_type value_type;
@@ -281,7 +281,7 @@
 // template function to call gees, default workspace type
 template< typename MatrixA, typename VectorWR, typename VectorWI,
         typename MatrixVS >
-inline integer_t gees( const char jobvs, const char sort,
+inline integer_t gees_2( const char jobvs, const char sort,
         logical_t* select, MatrixA& a, integer_t& sdim, VectorWR& wr,
         VectorWI& wi, MatrixVS& vs ) {
     typedef typename traits::matrix_traits< MatrixA >::value_type value_type;

Modified: sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_gees.cpp
==============================================================================
--- sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_gees.cpp	(original)
+++ sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_gees.cpp	2009-12-21 08:47:04 EST (Mon, 21 Dec 2009)
@@ -8,10 +8,13 @@
 
 #include "../../blas/test/random.hpp"
 
-#include <boost/numeric/bindings/lapack/gees.hpp>
+#include <boost/numeric/bindings/lapack/driver/gees.hpp>
 #include <boost/numeric/bindings/traits/ublas_matrix.hpp>
 #include <boost/numeric/bindings/traits/ublas_vector.hpp>
+#include <boost/numeric/bindings/traits/detail/utils.hpp>
 #include <boost/numeric/ublas/io.hpp>
+#include <boost/type_traits/is_complex.hpp>
+#include <boost/mpl/if.hpp>
 
 #include <iostream>
 #include <limits>
@@ -19,6 +22,34 @@
 
 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 MatrixVS,
+        typename Workspace >
+  static inline integer_t gees( const char jobvs, const char sort,
+        logical_t* select, MatrixA& a, integer_t& sdim, VectorW& w,
+        MatrixVS& vs, Workspace work ) {
+    typedef typename traits::matrix_traits< MatrixA >::value_type value_type;
+    traits::detail::array<value_type> wr(traits::vector_size(w));
+    traits::detail::array<value_type> wi(traits::vector_size(w));
+    return lapack::gees_2( jobvs, sort, select, a, sdim, wr, wi, vs, work );
+    traits::detail::interlace(traits::vector_storage(wr),
+                              traits::vector_storage(wr)+traits::vector_size(w),
+                              traits::vector_storage(wi),
+                              traits::vector_storage(w));
+  }
+};
+
+struct apply_complex {
+  template< typename MatrixA, typename VectorW, typename MatrixVS,
+        typename Workspace >
+  static inline integer_t gees( const char jobvs, const char sort,
+        logical_t* select, MatrixA& a, integer_t& sdim, VectorW& w,
+        MatrixVS& vs, Workspace work ) {
+    return lapack::gees( jobvs, sort, select, a, sdim, w, vs, work );
+  }
+};
 
 
 // Randomize a matrix
@@ -42,6 +73,7 @@
 
 template <typename T, typename W>
 int do_memory_type(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 std::complex< real_type >                                            complex_type ;
 
@@ -57,15 +89,17 @@
 
    randomize( a );
    matrix_type a2( a );
-
+   logical_t* select = 0;
+   integer_t sdim_info(0);
    // Compute Schur decomposition.
-   lapack::gees( a, e1, z, workspace ) ;
+   apply_t::gees( 'V', 'N', select, a, sdim_info, e1, z, workspace ) ;
 
    // Check Schur factorization
    if (norm_frobenius( prod( a2, z ) - prod( z, a ) )
            >= safety_factor*10.0* norm_frobenius( a2 ) * std::numeric_limits< real_type >::epsilon() ) return 255 ;
 
-   lapack::gees( a2, e2, workspace ) ;
+   matrix_type z_dummy( 1, 1 );
+   apply_t::gees( 'N', 'N', select, a2, sdim_info, e2, z_dummy, workspace ) ;
    if (norm_2( e1 - e2 ) > safety_factor*norm_2( e1 ) * std::numeric_limits< real_type >::epsilon()) return 255 ;
 
    if (norm_frobenius( a2 - a )
@@ -79,17 +113,20 @@
 template <typename T>
 struct Workspace {
    typedef ublas::vector<T>                         array_type ;
-   typedef lapack::detail::workspace1< array_type > type ;
+   typedef ublas::vector< bool >                    bool_array_type ;
+   typedef lapack::detail::workspace2< array_type,bool_array_type > type ;
 
    Workspace(size_t n)
    : work_( 3*n )
+   , bwork_(1)
    {}
 
    type operator() () {
-      return lapack::workspace(work_) ;
+      return lapack::workspace(work_, bwork_) ;
    }
 
    array_type work_ ;
+   bool_array_type bwork_;
 };
 
 
@@ -97,19 +134,22 @@
 struct Workspace< std::complex<T> > {
    typedef ublas::vector<T>                                                 real_array_type ;
    typedef ublas::vector< std::complex<T> >                                 complex_array_type ;
-   typedef lapack::detail::workspace2< complex_array_type,real_array_type > type ;
+   typedef ublas::vector< bool >                                            bool_array_type ;
+   typedef lapack::detail::workspace3< complex_array_type,real_array_type,bool_array_type > type ;
 
    Workspace(size_t n)
    : work_( 2*n )
    , rwork_( n )
+   , bwork_(1)
    {}
 
    type operator() () {
-      return lapack::workspace(work_, rwork_) ;
+      return lapack::workspace(work_, rwork_, bwork_) ;
    }
 
    complex_array_type work_ ;
    real_array_type    rwork_ ;
+   bool_array_type    bwork_;
 };
 
 
@@ -121,7 +161,7 @@
    if (do_memory_type<T,lapack::minimal_workspace>( n, lapack::minimal_workspace() ) ) return 255 ;
 
    Workspace<T> work( n );
-   do_memory_type<T,typename Workspace<T>::type >( n, work() );
+   if (do_memory_type<T,typename Workspace<T>::type >( n, work() ) ) return 255 ;
    return 0;
 } // do_value_type()
 

Modified: sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_geev.cpp
==============================================================================
--- sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_geev.cpp	(original)
+++ sandbox/numeric_bindings/libs/numeric/bindings/lapack/test/ublas_geev.cpp	2009-12-21 08:47:04 EST (Mon, 21 Dec 2009)
@@ -6,7 +6,7 @@
 #include <boost/numeric/ublas/io.hpp>
 #include <boost/numeric/bindings/traits/ublas_matrix.hpp>
 #include <boost/numeric/bindings/traits/ublas_vector.hpp>
-#include <boost/numeric/bindings/lapack/geev.hpp>
+#include <boost/numeric/bindings/lapack/driver/geev.hpp>
 #include "utils.h"
 
 using std::cout;
@@ -30,15 +30,18 @@
 }
 void geev(int n){
     cout << "\nCalculating eigenvalues using LAPACK's geev." << endl;
-    ublas::matrix<double, ublas::column_major> A(n,n);
+    //ublas::matrix<double, ublas::column_major> A(n,n);
+    //falling back to fully complex case for the moment,
+    //because the reassembly of the result in the real case is a bit complicated.
+    ublas::matrix<complex<double>, ublas::column_major> A(n,n);
     Hessenberg(A);
     print_m(A);
 
     ublas::vector<complex<double> > values(n);
-    ublas::matrix<complex<double>, ublas::column_major>* Vectors_left = 0;
+    ublas::matrix<complex<double>, ublas::column_major> Vectors_left(1,1);
     ublas::matrix<complex<double>, ublas::column_major> Vectors_right(n,n);
 
-    lapack::geev(A, values, Vectors_left, &Vectors_right, lapack::optimal_workspace());
+    lapack::geev('N','V', A, values, Vectors_left, Vectors_right, lapack::optimal_workspace());
     print_v(values, "values"); cout << endl;
     print_m(Vectors_right, "Vectors_right"); cout << endl;