[Boost-commit] svn:boost r58552 - sandbox/numeric_bindings/libs/numeric/bindings/blas/test

Subject: [Boost-commit] svn:boost r58552 - sandbox/numeric_bindings/libs/numeric/bindings/blas/test
From: thomas.klimpel_at_[hidden]
Date: 2009-12-28 18:52:34

Author: klimpel
Date: 2009-12-28 18:52:34 EST (Mon, 28 Dec 2009)
New Revision: 58552
URL: http://svn.boost.org/trac/boost/changeset/58552

Log:
Continue of regression tests update
- complex versions of asum and nrm2 are missing, use manual emulation as temporary substitute in blas1.cpp
Text files modified: 
   sandbox/numeric_bindings/libs/numeric/bindings/blas/test/blas1.cpp |    53 +++++++++++++++++++++++++++++-          
   sandbox/numeric_bindings/libs/numeric/bindings/blas/test/blas2.cpp |     3 +                                       
   sandbox/numeric_bindings/libs/numeric/bindings/blas/test/blas3.cpp |    67 ++++++++++++++++++++++++++++++++++----- 
   3 files changed, 109 insertions(+), 14 deletions(-)

Modified: sandbox/numeric_bindings/libs/numeric/bindings/blas/test/blas1.cpp
==============================================================================
--- sandbox/numeric_bindings/libs/numeric/bindings/blas/test/blas1.cpp	(original)
+++ sandbox/numeric_bindings/libs/numeric/bindings/blas/test/blas1.cpp	2009-12-28 18:52:34 EST (Mon, 28 Dec 2009)
@@ -10,10 +10,20 @@
 
 #include <boost/numeric/ublas/vector_proxy.hpp>
 #include <boost/numeric/ublas/matrix_proxy.hpp>
+#include <boost/numeric/ublas/blas.hpp>
 #include <boost/numeric/bindings/traits/ublas_vector.hpp>
 #include <boost/numeric/bindings/traits/ublas_matrix.hpp>
 #include <boost/numeric/bindings/traits/std_vector.hpp>
-#include <boost/numeric/bindings/blas/blas1.hpp>
+#include <boost/numeric/bindings/blas/level1/axpy.hpp>
+#include <boost/numeric/bindings/blas/level1/asum.hpp>
+#include <boost/numeric/bindings/blas/level1/copy.hpp>
+#include <boost/numeric/bindings/blas/level1/dot.hpp>
+#include <boost/numeric/bindings/blas/level1/dotc.hpp>
+#include <boost/numeric/bindings/blas/level1/dotu.hpp>
+#include <boost/numeric/bindings/blas/level1/nrm2.hpp>
+#include <boost/numeric/bindings/blas/level1/scal.hpp>
+#include <boost/type_traits/is_complex.hpp>
+#include <boost/mpl/if.hpp>
 
 #include <vector>
 #include <complex>
@@ -21,6 +31,42 @@
 #include <limits>
 #include <cmath>
 
+namespace traits = boost::numeric::bindings::traits;
+
+struct apply_real {
+  template< typename VectorX >
+  static inline typename traits::type_traits< typename traits::vector_traits< VectorX >::value_type >::real_type
+  asum( const VectorX& x ) {
+    return boost::numeric::bindings::blas::asum( x );
+  }
+  template< typename VectorX >
+  static inline typename traits::type_traits< typename traits::vector_traits< VectorX >::value_type >::real_type
+  nrm2( const VectorX& x ) {
+    return boost::numeric::bindings::blas::nrm2( x );
+  }
+};
+
+struct apply_complex {
+  template< typename VectorX >
+  static inline typename traits::type_traits< typename traits::vector_traits< VectorX >::value_type >::real_type
+  asum( const VectorX& x ) {
+    return abs_sum( x );
+  }
+  template< typename VectorX >
+  static inline typename traits::type_traits< typename traits::vector_traits< VectorX >::value_type >::real_type
+  nrm2( const VectorX& x ) {
+    using namespace std;
+    typedef typename traits::vector_traits< VectorX >::value_type value_type;
+    typedef typename traits::type_traits< value_type >::real_type real_type;
+    real_type t = real_type();
+    for (size_t i = 0, s = x.size(); i < s; ++ i) {
+      real_type u = norm(x[i]);
+      t += u;
+    }
+    return sqrt (t);
+  }
+};
+
 
 
 // Randomize a vector (using functions from random.hpp)
@@ -77,6 +123,7 @@
 
   template <typename V>
   int operator()(V& v) const {
+     typedef typename boost::mpl::if_<boost::is_complex<T>, apply_complex, apply_real>::type apply_t;
      using namespace boost::numeric::bindings::blas ;
 
      typedef typename V::value_type                                                        value_type ;
@@ -87,13 +134,13 @@
         v[i] = v_ref_(i);
 
      // Test blas routines and compare with reference
-     real_type nrm = nrm2( v );
+     real_type nrm = apply_t::nrm2( v );
      if ( std::abs(nrm - norm_2(v_ref_)) > std::numeric_limits< real_type >::epsilon() * norm_2(v_ref_)) {
        std::cout << "nrm2 : " << std::abs(nrm - norm_2(v_ref_)) << " > " << std::numeric_limits< real_type >::epsilon() * norm_2(v_ref_) << std::endl ;
        return 255 ;
      }
 
-     nrm = asum( v );
+     nrm = apply_t::asum( v );
      if ( std::abs(nrm - abs_sum(v_ref_)) > std::numeric_limits< real_type >::epsilon() * abs_sum(v_ref_)) {
        std::cout << "asum : " << std::abs(nrm - abs_sum(v_ref_)) << " > " << std::numeric_limits< real_type >::epsilon() * abs_sum(v_ref_) << std::endl ;
        return 255 ;

Modified: sandbox/numeric_bindings/libs/numeric/bindings/blas/test/blas2.cpp
==============================================================================
--- sandbox/numeric_bindings/libs/numeric/bindings/blas/test/blas2.cpp	(original)
+++ sandbox/numeric_bindings/libs/numeric/bindings/blas/test/blas2.cpp	2009-12-28 18:52:34 EST (Mon, 28 Dec 2009)
@@ -13,7 +13,8 @@
 #include <iomanip>
 
 #include "blas.hpp"
-#include <boost/numeric/bindings/blas/blas2.hpp>
+#include <boost/numeric/bindings/blas/level2/gemv.hpp>
+#include <boost/numeric/bindings/traits/transpose.hpp>
 #include <boost/numeric/bindings/traits/ublas_vector.hpp>
 #include <boost/numeric/bindings/traits/ublas_matrix.hpp>
 

Modified: sandbox/numeric_bindings/libs/numeric/bindings/blas/test/blas3.cpp
==============================================================================
--- sandbox/numeric_bindings/libs/numeric/bindings/blas/test/blas3.cpp	(original)
+++ sandbox/numeric_bindings/libs/numeric/bindings/blas/test/blas3.cpp	2009-12-28 18:52:34 EST (Mon, 28 Dec 2009)
@@ -13,17 +13,49 @@
 #include <boost/numeric/bindings/traits/ublas_matrix.hpp>
 #include <boost/numeric/bindings/traits/ublas_symmetric.hpp>
 #include <boost/numeric/bindings/traits/ublas_hermitian.hpp>
-#include <boost/numeric/bindings/blas/blas3.hpp>
+#include <boost/numeric/bindings/blas/level3/gemm.hpp>
+#include <boost/numeric/bindings/blas/level3/herk.hpp>
+#include <boost/numeric/bindings/blas/level3/syrk.hpp>
+#include <boost/numeric/bindings/blas/level3/trsm.hpp>
 #include <boost/numeric/ublas/matrix.hpp>
 #include <boost/numeric/ublas/matrix_proxy.hpp>
 #include <boost/numeric/ublas/symmetric.hpp>
 #include <boost/numeric/ublas/hermitian.hpp>
 #include <boost/numeric/ublas/triangular.hpp>
 #include <boost/numeric/ublas/io.hpp>
+#include <boost/type_traits/is_complex.hpp>
+#include <boost/mpl/if.hpp>
 
 #include <iostream>
 #include <complex>
 
+namespace traits = boost::numeric::bindings::traits;
+
+struct apply_real {
+  template< typename MatrixA, typename MatrixC >
+  static inline typename boost::numeric::bindings::blas::herk_impl< typename traits::matrix_traits<
+        MatrixA >::value_type >::return_type
+  herk( const char trans, const typename traits::type_traits<
+        typename traits::matrix_traits<
+        MatrixA >::value_type >::real_type alpha, const MatrixA& a,
+        const typename traits::type_traits< typename traits::matrix_traits<
+        MatrixA >::value_type >::real_type beta, MatrixC& c ) {
+    return boost::numeric::bindings::blas::syrk( trans, alpha, a, beta, c );
+  }
+};
+
+struct apply_complex {
+  template< typename MatrixA, typename MatrixC >
+  static inline typename boost::numeric::bindings::blas::herk_impl< typename traits::matrix_traits<
+        MatrixA >::value_type >::return_type
+  herk( const char trans, const typename traits::type_traits<
+        typename traits::matrix_traits<
+        MatrixA >::value_type >::real_type alpha, const MatrixA& a,
+        const typename traits::type_traits< typename traits::matrix_traits<
+        MatrixA >::value_type >::real_type beta, MatrixC& c ) {
+    return boost::numeric::bindings::blas::herk( trans, alpha, a, beta, c );
+  }
+};
 
 // Randomize a matrix
 template <typename M>
@@ -105,6 +137,7 @@
    typedef boost::numeric::ublas::matrix<value_type, boost::numeric::ublas::column_major>  ref_matrix_type ;
    typedef typename boost::numeric::bindings::traits::type_traits< value_type >::real_type real_type ;
    typedef typename ref_matrix_type::size_type                                             size_type ;
+   typedef typename boost::mpl::if_<boost::is_complex<value_type>, apply_complex, apply_real>::type apply_t;
 
    Syrk2(const M1& a,
         const ref_matrix_type& a_ref)
@@ -123,24 +156,30 @@
       real_type beta = -3.0 ;
 
       c.assign( c_ref_ );
-      boost::numeric::bindings::blas::syrk( 'U', 'N', value_type(alpha), a_, value_type(beta), c ) ;
+      //boost::numeric::bindings::blas::syrk( 'U', 'N', value_type(alpha), a_, value_type(beta), c ) ;
+      symmetric_adaptor<M, upper> sc( c );
+      boost::numeric::bindings::blas::syrk( 'N', value_type(alpha), a_, value_type(beta), sc ) ;
       if ( norm_frobenius( upper_part( c - (beta*c_ref_ + alpha * prod( a_ref_, trans( a_ref_ ) ) ) ) )
           > std::numeric_limits< real_type >::epsilon() * norm_frobenius( upper_part(c) ) ) return 255;
 
       c.assign( c_ref_ );
       symmetric_adaptor<M, upper> c_s( c );
-      boost::numeric::bindings::blas::syrk( 'U', 'N', value_type(alpha), a_, value_type(beta), c_s ) ;
+      //boost::numeric::bindings::blas::syrk( 'U', 'N', value_type(alpha), a_, value_type(beta), c_s ) ;
+      boost::numeric::bindings::blas::syrk( 'N', value_type(alpha), a_, value_type(beta), c_s ) ;
       if ( norm_frobenius( upper_part( c_s - (beta*c_ref_ + alpha * prod( a_ref_, trans( a_ref_ ) ) ) ) )
           > std::numeric_limits< real_type >::epsilon() * norm_frobenius( upper_part(c_s) ) ) return 255;
 
       c.assign( c_ref_ );
-      boost::numeric::bindings::blas::herk( 'U', 'N', alpha, a_, beta, c ) ;
+      //boost::numeric::bindings::blas::herk( 'U', 'N', alpha, a_, beta, c ) ;
+      hermitian_adaptor<M, upper> hc( c );
+      apply_t::herk( 'N', alpha, a_, beta, hc ) ;
       if ( norm_frobenius( upper_part( c - (beta*c_ref_ + alpha * prod( a_ref_, herm( a_ref_ ) ) ) ) )
           > std::numeric_limits< real_type >::epsilon() * norm_frobenius( upper_part(c) ) ) return 255;
 
       c.assign( c_ref_ );
       hermitian_adaptor<M, upper> c_h( c );
-      boost::numeric::bindings::blas::herk( 'U', 'N', alpha, a_, beta, c_h ) ;
+      //boost::numeric::bindings::blas::herk( 'U', 'N', alpha, a_, beta, c_h ) ;
+      apply_t::herk( 'N', alpha, a_, beta, c_h ) ;
       if ( norm_frobenius( upper_part( c_h - (beta*c_ref_ + alpha * prod( a_ref_, herm( a_ref_ ) ) ) ) )
           > std::numeric_limits< real_type >::epsilon() * norm_frobenius( upper_part(c_h) ) ) return 255;
 
@@ -304,27 +343,35 @@
       real_type alpha = 2.0 ;
 
       a.assign( a_ref_ );
-      boost::numeric::bindings::blas::trsm( 'L', 'U', 'N', 'N', value_type(alpha), b_ref_, a ) ;
+      //boost::numeric::bindings::blas::trsm( 'L', 'U', 'N', 'N', value_type(alpha), b_ref_, a ) ;
+      symmetric_adaptor<const ref_matrix_type, upper> sb_ref( b_ref_ );
+      boost::numeric::bindings::blas::trsm( 'L', 'N', 'N', value_type(alpha), sb_ref, a ) ;
       if ( norm_frobenius( alpha * a_ref_ - prod( upper_part(b_ref_), a ) )
           > std::numeric_limits< real_type >::epsilon() * norm_frobenius( upper_part(b_ref_) ) * norm_frobenius( a) ) return 255;
 
       a.assign( a_ref_ );
-      boost::numeric::bindings::blas::trsm( 'R', 'U', 'N', 'N', value_type(alpha), c_, a ) ;
+      //boost::numeric::bindings::blas::trsm( 'R', 'U', 'N', 'N', value_type(alpha), c_, a ) ;
+      symmetric_adaptor<const M1, upper> suc( c_ );
+      boost::numeric::bindings::blas::trsm( 'R', 'N', 'N', value_type(alpha), suc, a ) ;
       if ( norm_frobenius( alpha * a_ref_ - prod( a, upper_part(c_) ) )
           > std::numeric_limits< real_type >::epsilon() * norm_frobenius( upper_part(c_) ) * norm_frobenius( a ) ) return 255;
 
       a.assign( a_ref_ );
-      boost::numeric::bindings::blas::trsm( 'R', 'L', 'N', 'N', value_type(alpha), c_, a ) ;
+      //boost::numeric::bindings::blas::trsm( 'R', 'L', 'N', 'N', value_type(alpha), c_, a ) ;
+      symmetric_adaptor<const M1, lower> slc( c_ );
+      boost::numeric::bindings::blas::trsm( 'R', 'N', 'N', value_type(alpha), slc, a ) ;
       if ( norm_frobenius( alpha * a_ref_ - prod( a, lower_part(c_) ) )
           > std::numeric_limits< real_type >::epsilon() * norm_frobenius( lower_part(c_) ) * norm_frobenius( a ) ) return 255;
 
       a.assign( a_ref_ );
-      boost::numeric::bindings::blas::trsm( 'L', 'U', 'T', 'N', value_type(alpha), b_ref_, a ) ;
+      //boost::numeric::bindings::blas::trsm( 'L', 'U', 'T', 'N', value_type(alpha), b_ref_, a ) ;
+      boost::numeric::bindings::blas::trsm( 'L', 'T', 'N', value_type(alpha), sb_ref, a ) ;
       if ( norm_frobenius( alpha * a_ref_ - prod( trans(upper_part(b_ref_)), a ) )
           > std::numeric_limits< real_type >::epsilon() * norm_frobenius( upper_part(b_ref_) ) * norm_frobenius( a) ) return 255;
 
       a.assign( a_ref_ );
-      boost::numeric::bindings::blas::trsm( 'L', 'U', 'N', 'U', value_type(alpha), b_ref_, a ) ;
+      //boost::numeric::bindings::blas::trsm( 'L', 'U', 'N', 'U', value_type(alpha), b_ref_, a ) ;
+      boost::numeric::bindings::blas::trsm( 'L', 'N', 'U', value_type(alpha), sb_ref, a ) ;
       if ( norm_frobenius( alpha * a_ref_ - prod( unit_upper_part(b_ref_), a ) )
           > std::numeric_limits< real_type >::epsilon() * norm_frobenius( unit_upper_part(b_ref_) ) * norm_frobenius( a) ) return 255;