blas1.hpp

Go to the documentation of this file.
00001 //
00002 //  Copyright Toon Knapen and Kresimir Fresl
00003 //
00004 // Distributed under the Boost Software License, Version 1.0.
00005 // (See accompanying file LICENSE_1_0.txt or copy at
00006 // http://www.boost.org/LICENSE_1_0.txt)
00007 //
00008 
00009 #ifndef BOOST_NUMERIC_BINDINGS_BLAS_BLAS1_HPP
00010 #define BOOST_NUMERIC_BINDINGS_BLAS_BLAS1_HPP
00011 
00012 #include <boost/numeric/bindings/blas/blas1_overloads.hpp>
00013 #include <boost/numeric/bindings/traits/vector_traits.hpp>
00014 #include <boost/static_assert.hpp>
00015 #include <boost/type_traits/is_same.hpp>
00016 #include <cassert> 
00017 
00018 namespace boost { namespace numeric { namespace bindings { namespace blas {
00019 
00020   // x <- y
00021   template < typename vector_x_type, typename vector_y_type >
00022   void copy(const vector_x_type &x, vector_y_type &y )
00023   {
00024 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00025     BOOST_STATIC_ASSERT( ( boost::is_same< typename traits::vector_traits<vector_x_type>::value_type, typename traits::vector_traits<vector_y_type>::value_type >::value ) ) ;
00026 #else
00027     BOOST_STATIC_ASSERT( ( boost::is_same< typename vector_x_type::value_type, typename vector_y_type::value_type >::value ) ) ;
00028 #endif
00029 
00030     const int n =  traits::vector_size( x ) ;
00031     assert( n==traits::vector_size( y ) ) ;
00032     const int stride_x = traits::vector_stride( x ) ;
00033     const int stride_y = traits::vector_stride( y ) ;
00034     typename traits::vector_traits<vector_x_type>::value_type const *x_ptr = traits::vector_storage( x ) ;
00035     typename traits::vector_traits<vector_y_type>::value_type *y_ptr = traits::vector_storage( y ) ;
00036 
00037     detail::copy( n, x_ptr, stride_x, y_ptr, stride_y ) ;
00038   }
00039 
00040 
00041   // x <- alpha * x
00042   template < typename value_type, typename vector_type >
00043   void scal(const value_type &alpha, vector_type &x )
00044   {
00045 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00046     BOOST_STATIC_ASSERT( ( boost::is_same< value_type, typename traits::vector_traits<vector_type>::value_type >::value ) ) ;
00047 #else
00048     BOOST_STATIC_ASSERT( ( boost::is_same< value_type, typename vector_type::value_type >::value ) ) ;
00049 #endif
00050 
00051     const int n =  traits::vector_size( x ) ;
00052     const int stride = traits::vector_stride( x ) ;
00053     value_type *x_ptr = traits::vector_storage( x ) ;
00054 
00055     detail::scal( n, alpha, x_ptr, stride ) ;
00056   }
00057 
00058 
00059   // y <- alpha * x + y
00060   template < typename value_type, typename vector_type_x, typename vector_type_y >
00061   void axpy(const value_type& alpha, const vector_type_x &x, vector_type_y &y )
00062   { 
00063 #ifdef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00064     BOOST_STATIC_ASSERT( ( is_same< value_type, typename vector_type_x::value_type >::value ) ) ;
00065     BOOST_STATIC_ASSERT( ( is_same< value_type, typename vector_type_y::value_type >::value ) ) ;
00066 #else
00067     BOOST_STATIC_ASSERT( ( is_same< value_type, typename traits::vector_traits< vector_type_x >::value_type >::value ) ) ;
00068     BOOST_STATIC_ASSERT( ( is_same< value_type, typename traits::vector_traits< vector_type_y >::value_type >::value ) ) ;
00069 #endif
00070     assert( traits::vector_size( x ) == traits::vector_size( y ) ) ;
00071 
00072     const int n = traits::vector_size( x ) ;
00073     const int stride_x = traits::vector_stride( x ) ;
00074     const int stride_y = traits::vector_stride( y ) ;
00075     const value_type *x_ptr = traits::vector_storage( x ) ;
00076     value_type *y_ptr = traits::vector_storage( y ) ;
00077 
00078     detail::axpy( n, alpha, x_ptr, stride_x, y_ptr, stride_y ) ; 
00079   }
00080 
00081 
00082   // dot <- x^T * y  (real vectors)
00083   template < typename vector_type_x, typename vector_type_y >
00084 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00085   typename traits::vector_traits< vector_type_x >::value_type 
00086 #else
00087   typename vector_type_x::value_type
00088 #endif
00089   dot(const vector_type_x &x, const vector_type_y &y)
00090   {
00091 #ifdef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00092     BOOST_STATIC_ASSERT( ( is_same< typename vector_type_y::value_type, typename vector_type_x::value_type >::value ) ) ;
00093 #else
00094     BOOST_STATIC_ASSERT( ( is_same< typename traits::vector_traits< vector_type_y >::value_type, typename traits::vector_traits< vector_type_x >::value_type >::value ) ) ;
00095 #endif
00096 
00097     assert( traits::vector_size( x ) == traits::vector_size( y ) ) ;
00098 
00099     typedef
00100 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00101     typename traits::vector_traits< vector_type_x >::value_type 
00102 #else
00103     typename vector_type_x::value_type
00104 #endif
00105     value_type ;
00106 
00107     const int n = traits::vector_size( x ) ;
00108     const int stride_x = traits::vector_stride( x ) ;
00109     const int stride_y = traits::vector_stride( y ) ;
00110 
00111     const value_type *x_ptr = traits::vector_storage( x ) ;
00112     const value_type *y_ptr = traits::vector_storage( y ) ;
00113 
00114     return detail::dot( n, x_ptr, stride_x, y_ptr, stride_y ) ;
00115   }
00116 
00117   // dotu <- x^T * y  (complex vectors)
00118   template < typename vector_type_x, typename vector_type_y >
00119 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00120   typename traits::vector_traits< vector_type_x >::value_type 
00121 #else
00122   typename vector_type_x::value_type
00123 #endif
00124   dotu(const vector_type_x &x, const vector_type_y &y)
00125   {
00126 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00127     BOOST_STATIC_ASSERT( ( is_same< typename traits::vector_traits< vector_type_y >::value_type, typename traits::vector_traits< vector_type_x >::value_type >::value ) ) ;
00128 #else
00129     BOOST_STATIC_ASSERT( ( is_same< typename vector_type_y::value_type, typename vector_type_x::value_type >::value ) ) ;
00130 #endif
00131     assert( traits::vector_size( x ) == traits::vector_size( y ) ) ;
00132 
00133     typedef typename vector_type_x::value_type value_type ;
00134 
00135     const int n = traits::vector_size( x ) ;
00136     const int stride_x = traits::vector_stride( x ) ;
00137     const int stride_y = traits::vector_stride( y ) ;
00138     const value_type *x_ptr = traits::vector_storage( x ) ;
00139     const value_type *y_ptr = traits::vector_storage( y ) ;
00140     
00141     value_type ret ;
00142     detail::dotu( ret, n, x_ptr, stride_x, y_ptr, stride_y ) ;
00143     return ret;
00144   }
00145 
00146   // dotc <- x^H * y  (complex vectors) 
00147   template < typename vector_type_x, typename vector_type_y >
00148 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00149   typename traits::vector_traits< vector_type_x >::value_type 
00150 #else
00151   typename vector_type_x::value_type
00152 #endif
00153   dotc(const vector_type_x &x, const vector_type_y &y)
00154   {
00155 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00156     BOOST_STATIC_ASSERT( ( is_same< typename traits::vector_traits< vector_type_y >::value_type, typename traits::vector_traits< vector_type_x >::value_type >::value ) ) ;
00157 #else
00158     BOOST_STATIC_ASSERT( ( is_same< typename vector_type_y::value_type, typename vector_type_x::value_type >::value ) ) ;
00159 #endif
00160     assert( traits::vector_size( x ) == traits::vector_size( y ) ) ;
00161 
00162     typedef typename vector_type_x::value_type value_type ;
00163 
00164     const int n = traits::vector_size( x ) ;
00165     const int stride_x = traits::vector_stride( x ) ;
00166     const int stride_y = traits::vector_stride( y ) ;
00167     const value_type *x_ptr = traits::vector_storage( x ) ;
00168     const value_type *y_ptr = traits::vector_storage( y ) ;
00169     
00170     value_type ret ;
00171     detail::dotc( ret, n, x_ptr, stride_x, y_ptr, stride_y ) ;
00172     return ret;
00173   }
00174 
00175 
00176   // nrm2 <- ||x||_2
00177   template < typename vector_type >
00178 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00179   typename traits::type_traits< typename traits::vector_traits< vector_type >::value_type >::real_type
00180 #else
00181   typename traits::type_traits< typename vector_type::value_type >::real_type
00182 #endif
00183   nrm2(const vector_type &x) 
00184   {
00185 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00186     typedef typename traits::vector_traits< vector_type >::value_type value_type;
00187 #else
00188     typedef vector_type::value_type value_type ;
00189 #endif
00190     const int n = traits::vector_size( x ) ;
00191     const int stride_x = traits::vector_stride( x ) ;
00192     const value_type *x_ptr = traits::vector_storage( x ) ;
00193 
00194     return detail::nrm2( n, x_ptr, stride_x ) ;
00195   }
00196 
00197 
00198 
00199   // asum <- ||x||_1
00200   // .. for now works only with real vectors
00201   template < typename vector_type >
00202 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00203   typename traits::type_traits< typename traits::vector_traits< vector_type >::value_type >::real_type 
00204 #else
00205   typename traits::type_traits< typename vector_type::value_type >::real_type
00206 #endif
00207   asum(const vector_type &x) 
00208   {
00209 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00210     typedef typename traits::vector_traits< vector_type >::value_type value_type;
00211 #else
00212     typedef vector_type::value_type value_type ;
00213 #endif
00214 
00215     const int n = traits::vector_size( x ) ;
00216     const int stride_x = traits::vector_stride( x ) ;
00217     const value_type *x_ptr = traits::vector_storage( x ) ;
00218 
00219     return detail::asum( n, x_ptr, stride_x ) ;
00220   }
00221 
00222 }}}}
00223 
00224 #endif // BOOST_NUMERIC_BINDINGS_BLAS_BLAS1_HPP