cblas1.hpp

Go to the documentation of this file.
00001 /*
00002  * 
00003  * Copyright (c) Kresimir Fresl 2002 
00004  *
00005  * Distributed under the Boost Software License, Version 1.0.
00006  * (See accompanying file LICENSE_1_0.txt or copy at
00007  * http://www.boost.org/LICENSE_1_0.txt)
00008  *
00009  * Author acknowledges the support of the Faculty of Civil Engineering, 
00010  * University of Zagreb, Croatia.
00011  *
00012  */
00013 
00014 #ifndef BOOST_NUMERIC_BINDINGS_CBLAS_LEVEL_1_HPP
00015 #define BOOST_NUMERIC_BINDINGS_CBLAS_LEVEL_1_HPP
00016 
00017 #include <cassert>
00018 
00019 #include <boost/numeric/bindings/traits/type_traits.hpp>
00020 #include <boost/numeric/bindings/traits/vector_traits.hpp>
00021 #include <boost/numeric/bindings/atlas/cblas1_overloads.hpp>
00022 
00023 #ifndef BOOST_NUMERIC_BINDINGS_NO_TYPE_CHECK
00024 #  include <boost/type_traits/same_traits.hpp>
00025 #  include <boost/static_assert.hpp>
00026 #endif 
00027 
00028 namespace boost { namespace numeric { namespace bindings { 
00029 
00030   namespace atlas {
00031 
00032     // x_i <- alpha for all i
00033     template <typename T, typename Vct> 
00034     inline 
00035     void set (T const& alpha, Vct& x) {
00036       detail::set (traits::vector_size (x), alpha, 
00037                    traits::vector_storage (x), traits::vector_stride (x)); 
00038     }
00039 
00040     // y <- x
00041     template <typename VctX, typename VctY>
00042     inline 
00043     void copy (VctX const& x, VctY& y) {
00044       assert (traits::vector_size (y) >= traits::vector_size (x));
00045       detail::copy (traits::vector_size (x),
00046 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00047                     traits::vector_storage (x), 
00048 #else
00049                     traits::vector_storage_const (x), 
00050 #endif 
00051                     traits::vector_stride (x), 
00052                     traits::vector_storage (y), traits::vector_stride (y)); 
00053     }
00054 
00055     // x <-> y
00056     template <typename VctX, typename VctY>
00057     inline 
00058     void swap (VctX& x, VctY& y) {
00059       assert (traits::vector_size (y) >= traits::vector_size (x));
00060       detail::swap (traits::vector_size (x),
00061                     traits::vector_storage (x), traits::vector_stride (x), 
00062                     traits::vector_storage (y), traits::vector_stride (y)); 
00063     }
00064 
00065     // x <- alpha * x
00066     template <typename T, typename Vct> 
00067     inline 
00068     void scal (T const& alpha, Vct& x) {
00069 #ifndef BOOST_NUMERIC_BINDINGS_NO_TYPE_CHECK
00070       typedef traits::vector_traits<Vct> vtraits;
00071       BOOST_STATIC_ASSERT(
00072        (boost::is_same<T, typename vtraits::value_type>::value
00073         || 
00074         boost::is_same<T, 
00075           typename traits::type_traits<typename vtraits::value_type>::real_type
00076         >::value
00077         ));
00078 #endif 
00079       detail::scal (traits::vector_size (x), alpha, 
00080                     traits::vector_storage (x), traits::vector_stride (x)); 
00081     }
00082 
00083     // y <- alpha * x + y
00084     template <typename T, typename VctX, typename VctY>
00085     inline 
00086     void axpy (T const& alpha, VctX const& x, VctY& y) {
00087       assert (traits::vector_size (y) >= traits::vector_size (x));
00088       detail::axpy (traits::vector_size (x), alpha, 
00089 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00090                     traits::vector_storage (x), 
00091 #else
00092                     traits::vector_storage_const (x), 
00093 #endif
00094                     traits::vector_stride (x), 
00095                     traits::vector_storage (y), traits::vector_stride (y)); 
00096     }
00097 
00098     // y <- x + y
00099     template <typename VctX, typename VctY>
00100     inline 
00101     void xpy (VctX const& x, VctY& y) {
00102 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00103       typedef typename traits::vector_traits<VctX>::value_type val_t; 
00104 #else
00105       typedef typename VctX::value_type val_t; 
00106 #endif
00107       axpy ((val_t) 1, x, y); 
00108     }
00109 
00110     // y <- alpha * x + beta * y
00111     template <typename T, typename VctX, typename VctY>
00112     inline 
00113     void axpby (T const& alpha, VctX const& x, 
00114                 T const& beta, VctY& y) { 
00115       assert (traits::vector_size (y) >= traits::vector_size (x));
00116       detail::axpby (traits::vector_size (x), alpha, 
00117 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00118                      traits::vector_storage (x), 
00119 #else
00120                      traits::vector_storage_const (x), 
00121 #endif
00122                      traits::vector_stride (x), 
00123                      beta, 
00124                      traits::vector_storage (y), traits::vector_stride (y)); 
00125     }
00126 
00128 
00129     // dot <- x^T * y 
00130     // .. real & complex types
00131     template <typename VctX, typename VctY>
00132     inline 
00133 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00134     typename traits::vector_traits<VctX>::value_type 
00135 #else
00136     typename VctX::value_type 
00137 #endif
00138     dot (VctX const& x, VctY const& y) {
00139       assert (traits::vector_size (y) >= traits::vector_size (x));
00140       return detail::dot (traits::vector_size (x),
00141 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00142                           traits::vector_storage (x), 
00143 #else
00144                           traits::vector_storage_const (x), 
00145 #endif
00146                           traits::vector_stride (x), 
00147 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00148                           traits::vector_storage (y), 
00149 #else
00150                           traits::vector_storage_const (y), 
00151 #endif
00152                           traits::vector_stride (y)); 
00153     }
00154 
00155     // dot <- x^T * y 
00156     // .. float only -- with double accumulation
00157     template <typename VctX, typename VctY>
00158     inline 
00159     double dsdot (VctX const& x, VctY const& y) {
00160       assert (traits::vector_size (y) >= traits::vector_size (x));
00161       return cblas_dsdot (traits::vector_size (x),
00162 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00163                           traits::vector_storage (x), 
00164 #else
00165                           traits::vector_storage_const (x), 
00166 #endif
00167                           traits::vector_stride (x), 
00168 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00169                           traits::vector_storage (y), 
00170 #else
00171                           traits::vector_storage_const (y), 
00172 #endif
00173                           traits::vector_stride (y)); 
00174     }
00175 
00176     // apdot <- alpha + x^T * y    
00177     // .. float only -- computation uses double precision 
00178     template <typename VctX, typename VctY>
00179     inline 
00180     float sdsdot (float const alpha, VctX const& x, VctY const& y) {
00181       assert (traits::vector_size (y) >= traits::vector_size (x));
00182       return cblas_sdsdot (traits::vector_size (x), alpha, 
00183 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00184                            traits::vector_storage (x), 
00185 #else
00186                            traits::vector_storage_const (x), 
00187 #endif
00188                            traits::vector_stride (x), 
00189 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00190                            traits::vector_storage (y), 
00191 #else
00192                            traits::vector_storage_const (y), 
00193 #endif
00194                            traits::vector_stride (y)); 
00195     }
00196 
00197     // dotu <- x^T * y 
00198     // .. complex types only
00199     // .. function
00200     template <typename VctX, typename VctY>
00201     inline 
00202 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00203     typename traits::vector_traits<VctX>::value_type 
00204 #else
00205     typename VctX::value_type 
00206 #endif
00207     dotu (VctX const& x, VctY const& y) {
00208       assert (traits::vector_size (y) >= traits::vector_size (x));
00209 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00210       typename traits::vector_traits<VctX>::value_type val;
00211 #else
00212       typename VctX::value_type val; 
00213 #endif
00214       detail::dotu (traits::vector_size (x),
00215 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00216                     traits::vector_storage (x), 
00217 #else
00218                     traits::vector_storage_const (x), 
00219 #endif
00220                     traits::vector_stride (x), 
00221 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00222                     traits::vector_storage (y), 
00223 #else
00224                     traits::vector_storage_const (y), 
00225 #endif
00226                     traits::vector_stride (y), 
00227                     &val);
00228       return val; 
00229     }
00230     // .. procedure 
00231     template <typename VctX, typename VctY>
00232     inline 
00233     void dotu (VctX const& x, VctY const& y, 
00234 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00235                typename traits::vector_traits<VctX>::value_type& val
00236 #else
00237                typename VctX::value_type& val
00238 #endif
00239     ) {
00240       assert (traits::vector_size (y) >= traits::vector_size (x));
00241       detail::dotu (traits::vector_size (x),
00242 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00243                     traits::vector_storage (x), 
00244 #else
00245                     traits::vector_storage_const (x), 
00246 #endif
00247                     traits::vector_stride (x), 
00248 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00249                     traits::vector_storage (y), 
00250 #else
00251                     traits::vector_storage_const (y), 
00252 #endif
00253                     traits::vector_stride (y), 
00254                     &val);
00255     }
00256 
00257     // dotc <- x^H * y 
00258     // .. complex types only
00259     // .. function
00260     template <typename VctX, typename VctY>
00261     inline 
00262 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00263     typename traits::vector_traits<VctX>::value_type 
00264 #else
00265     typename VctX::value_type 
00266 #endif
00267     dotc (VctX const& x, VctY const& y) {
00268       assert (traits::vector_size (y) >= traits::vector_size (x));
00269 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00270       typename traits::vector_traits<VctX>::value_type val;
00271 #else
00272       typename VctX::value_type val; 
00273 #endif
00274       detail::dotc (traits::vector_size (x),
00275 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00276                     traits::vector_storage (x), 
00277 #else
00278                     traits::vector_storage_const (x), 
00279 #endif
00280                     traits::vector_stride (x), 
00281 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00282                     traits::vector_storage (y), 
00283 #else
00284                     traits::vector_storage_const (y), 
00285 #endif
00286                     traits::vector_stride (y),
00287                     &val);
00288       return val; 
00289     }
00290     // .. procedure 
00291     template <typename VctX, typename VctY>
00292     inline 
00293     void dotc (VctX const& x, VctY const& y, 
00294 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00295                typename traits::vector_traits<VctX>::value_type& val
00296 #else
00297                typename VctX::value_type& val
00298 #endif
00299     ) {
00300       assert (traits::vector_size (y) >= traits::vector_size (x));
00301       detail::dotc (traits::vector_size (x),
00302 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00303                     traits::vector_storage (x), 
00304 #else
00305                     traits::vector_storage_const (x), 
00306 #endif
00307                     traits::vector_stride (x), 
00308 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00309                     traits::vector_storage (y), 
00310 #else
00311                     traits::vector_storage_const (y), 
00312 #endif
00313                     traits::vector_stride (y), 
00314                     &val);
00315     }
00316 
00317     // nrm2 <- ||x||_2
00318     template <typename Vct> 
00319     inline 
00320     typename traits::type_traits<
00321 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00322       typename traits::vector_traits<Vct>::value_type 
00323 #else 
00324       typename Vct::value_type 
00325 #endif 
00326     >::real_type 
00327     nrm2 (Vct const& x) {
00328       return detail::nrm2 (traits::vector_size (x),
00329 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00330                            traits::vector_storage (x), 
00331 #else
00332                            traits::vector_storage_const (x), 
00333 #endif
00334                            traits::vector_stride (x)); 
00335     }
00336 
00337     // asum <- ||re (x)|| + ||im (x)||
00338     template <typename Vct> 
00339     inline 
00340     typename traits::type_traits<
00341 #ifndef BOOST_NUMERIC_BINDINGS_POOR_MANS_TRAITS
00342       typename traits::vector_traits<Vct>::value_type 
00343 #else 
00344       typename Vct::value_type 
00345 #endif 
00346     >::real_type  
00347     asum (Vct const& x) {
00348       return detail::asum (traits::vector_size (x),
00349 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00350                            traits::vector_storage (x), 
00351 #else
00352                            traits::vector_storage_const (x), 
00353 #endif
00354                            traits::vector_stride (x)); 
00355     }
00356 
00357     // iamax <- 1st i: max (|re (x_i)| + |im (x_i)|)
00358     template <typename Vct> 
00359     inline 
00360     CBLAS_INDEX iamax (Vct const& x) {
00361       return detail::iamax (traits::vector_size (x),
00362 #ifndef BOOST_NO_FUNCTION_TEMPLATE_ORDERING
00363                             traits::vector_storage (x), 
00364 #else
00365                             traits::vector_storage_const (x), 
00366 #endif
00367                             traits::vector_stride (x)); 
00368     }
00369 
00370     // TO DO: plane rotations 
00371 
00372   } // namespace atlas
00373 
00374 }}} 
00375 
00376 #endif // BOOST_NUMERIC_BINDINGS_CBLAS_LEVEL_1_HPP