external/mkl/mkl_operations.hpp

0001 /*
0002  [auto_generated]
0003  boost/numeric/odeint/external/mkl/mkl_operations.hpp
0004
0005  [begin_description]
0006  Wrapper classes for intel math kernel library types.
0007  Get a free, non-commercial download of MKL at
0008  http://software.intel.com/en-us/articles/non-commercial-software-download/
0009  [end_description]
0010
0011  Copyright 2010-2011 Mario Mulansky
0012  Copyright 2011-2013 Karsten Ahnert
0013
0014  Distributed under the Boost Software License, Version 1.0.
0015  (See accompanying file LICENSE_1_0.txt or
0016  copy at http://www.boost.org/LICENSE_1_0.txt)
0017  */
0018
0019
0020 #ifndef BOOST_NUMERIC_ODEINT_EXTERNAL_MKL_MKL_OPERATIONS_HPP_INCLUDED
0021 #define BOOST_NUMERIC_ODEINT_EXTERNAL_MKL_MKL_OPERATIONS_HPP_INCLUDED
0022
0023 #include <iostream>
0024
0025 #include <mkl_cblas.h>
0026 #include <boost/numeric/odeint/algebra/default_operations.hpp>
0027
0028 /* exemplary example for writing bindings to the Intel MKL library
0029  * see test/mkl for how to use mkl with odeint
0030  * this is a quick and dirty implementation showing the general possibility.
0031  * It works only with containers based on double and sequential memory allocation.
0032  */
0033
0034 namespace boost {
0035 namespace numeric {
0036 namespace odeint {
0037
0038 /* only defined for doubles */
0039 struct mkl_operations
0040 {
0041     //template< class Fac1 , class Fac2 > struct scale_sum2;
0042
0043     template< class F1 = double , class F2 = F1 >
0044     struct scale_sum2
0045     {
0046         typedef double Fac1;
0047         typedef double Fac2;
0048         const Fac1 m_alpha1;
0049         const Fac2 m_alpha2;
0050
0051         scale_sum2( const Fac1 alpha1 , const Fac2 alpha2 ) : m_alpha1( alpha1 ) , m_alpha2( alpha2 ) { }
0052
0053         template< class T1 , class T2 , class T3 >
0054         void operator()( T1 &t1 , const T2 &t2 , const T3 &t3) const
0055         {   // t1 = m_alpha1 * t2 + m_alpha2 * t3;
0056             // we get Containers that have size() and [i]-access
0057
0058             const int n = t1.size();
0059             //boost::numeric::odeint::copy( t1 , t3 );
0060             if( &(t2[0]) != &(t1[0]) )
0061             {
0062                 cblas_dcopy( n , &(t2[0]) , 1 , &(t1[0]) , 1 );
0063             }
0064             cblas_dscal( n , m_alpha1 , &(t1[0]) , 1 );
0065             cblas_daxpy( n , m_alpha2 , &(t3[0]) , 1 , &(t1[0]) , 1 );
0066             //daxpby( &n , &m_alpha2 , &(t3[0]) , &one , &m_alpha1 , &(t1[0]) , &one );
0067         }
0068     };
0069
0070     template< class F1 = double , class F2 = F1 , class F3 = F2 >
0071     struct scale_sum3
0072     {
0073         typedef double Fac1;
0074         typedef double Fac2;
0075         typedef double Fac3;
0076         const Fac1 m_alpha1;
0077         const Fac2 m_alpha2;
0078         const Fac3 m_alpha3;
0079
0080         scale_sum3( const Fac1 alpha1 , const Fac2 alpha2 , const Fac3 alpha3 )
0081             : m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) { }
0082
0083         template< class T1 , class T2 , class T3 , class T4 >
0084         void operator()( T1 &t1 , const T2 &t2 , const T3 &t3 , const T4 &t4 ) const
0085         {   // t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4;
0086             // we get Containers that have size() and [i]-access
0087
0088             const int n = t1.size();
0089             //boost::numeric::odeint::copy( t1 , t3 );
0090             if( &(t2[0]) != &(t1[0]) )
0091             {
0092                 cblas_dcopy( n , &(t2[0]) , 1 , &(t1[0]) , 1 );
0093             }
0094             cblas_dscal( n , m_alpha1 , &(t1[0]) , 1 );
0095             cblas_daxpy( n , m_alpha2 , &(t3[0]) , 1 , &(t1[0]) , 1 );
0096             //daxpby( &n , &m_alpha2 , &(t3[0]) , &one , &m_alpha1 , &(t1[0]) , &one );
0097             cblas_daxpy( n , m_alpha3 , &(t4[0]) , 1 , &(t1[0]) , 1 );
0098         }
0099     };
0100
0101     template< class F1 = double , class F2 = F1 , class F3 = F2 , class F4 = F3 >
0102     struct scale_sum4
0103     {
0104         typedef double Fac1;
0105         typedef double Fac2;
0106         typedef double Fac3;
0107         typedef double Fac4;
0108         const Fac1 m_alpha1;
0109         const Fac2 m_alpha2;
0110         const Fac3 m_alpha3;
0111         const Fac4 m_alpha4;
0112
0113         scale_sum4( const Fac1 alpha1 , const Fac2 alpha2 , const Fac3 alpha3 , const Fac4 alpha4 )
0114             : m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) , m_alpha4( alpha4 ) { }
0115
0116         template< class T1 , class T2 , class T3 , class T4 , class T5 >
0117         void operator()( T1 &t1 , const T2 &t2 , const T3 &t3 , const T4 &t4 , const T5 &t5 ) const
0118         {   // t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5;
0119             // we get Containers that have size() and [i]-access
0120
0121             const int n = t1.size();
0122             //boost::numeric::odeint::copy( t1 , t3 );
0123             if( &(t2[0]) != &(t1[0]) )
0124             {
0125                 cblas_dcopy( n , &(t2[0]) , 1 , &(t1[0]) , 1 );
0126             }
0127
0128             cblas_dscal( n , m_alpha1 , &(t1[0]) , 1 );
0129             cblas_daxpy( n , m_alpha2 , &(t3[0]) , 1 , &(t1[0]) , 1 );
0130             //daxpby( &n , &m_alpha2 , &(t3[0]) , &one , &m_alpha1 , &(t1[0]) , &one );
0131             cblas_daxpy( n , m_alpha3 , &(t4[0]) , 1 , &(t1[0]) , 1 );
0132             cblas_daxpy( n , m_alpha4 , &(t5[0]) , 1 , &(t1[0]) , 1 );
0133         }
0134     };
0135
0136
0137     template< class F1 = double , class F2 = F1 , class F3 = F2 , class F4 = F3 , class F5 = F4 >
0138     struct scale_sum5
0139     {
0140         typedef double Fac1;
0141         typedef double Fac2;
0142         typedef double Fac3;
0143         typedef double Fac4;
0144         typedef double Fac5;
0145         const Fac1 m_alpha1;
0146         const Fac2 m_alpha2;
0147         const Fac3 m_alpha3;
0148         const Fac4 m_alpha4;
0149         const Fac5 m_alpha5;
0150
0151         scale_sum5( const Fac1 alpha1 , const Fac2 alpha2 , const Fac3 alpha3 , const Fac4 alpha4 , const Fac5 alpha5 )
0152             : m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) , m_alpha4( alpha4 ) , m_alpha5( alpha5 ) { }
0153
0154         template< class T1 , class T2 , class T3 , class T4 , class T5 , class T6   >
0155         void operator()( T1 &t1 , const T2 &t2 , const T3 &t3 , const T4 &t4 , const T5 &t5 , const T6 &t6 ) const
0156         {   // t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5 + m_alpha5 * t6;
0157             // we get Containers that have size() and [i]-access
0158
0159             const int n = t1.size();
0160             //boost::numeric::odeint::copy( t1 , t3 );
0161             if( &(t2[0]) != &(t1[0]) )
0162             {
0163                 cblas_dcopy( n , &(t2[0]) , 1 , &(t1[0]) , 1 );
0164             }
0165
0166             cblas_dscal( n , m_alpha1 , &(t1[0]) , 1 );
0167             cblas_daxpy( n , m_alpha2 , &(t3[0]) , 1 , &(t1[0]) , 1 );
0168             //daxpby( &n , &m_alpha2 , &(t3[0]) , &one , &m_alpha1 , &(t1[0]) , &one );
0169             cblas_daxpy( n , m_alpha3 , &(t4[0]) , 1 , &(t1[0]) , 1 );
0170             cblas_daxpy( n , m_alpha4 , &(t5[0]) , 1 , &(t1[0]) , 1 );
0171             cblas_daxpy( n , m_alpha5 , &(t6[0]) , 1 , &(t1[0]) , 1 );
0172         }
0173     };
0174
0175 };
0176
0177 } // odeint
0178 } // numeric
0179 } // boost
0180
0181 #endif // BOOST_NUMERIC_ODEINT_EXTERNAL_MKL_MKL_OPERATIONS_HPP_INCLUDED