3 boost/numeric/odeint/external/mkl/mkl_operations.hpp
6 Wrapper classes for intel math kernel library types.
7 Get a free, non-commercial download of MKL at
8 http://software.intel.com/en-us/articles/non-commercial-software-download/
11 Copyright 2010-2011 Mario Mulansky
12 Copyright 2011-2013 Karsten Ahnert
14 Distributed under the Boost Software License, Version 1.0.
15 (See accompanying file LICENSE_1_0.txt or
16 copy at http://www.boost.org/LICENSE_1_0.txt)
20 #ifndef BOOST_NUMERIC_ODEINT_EXTERNAL_MKL_MKL_OPERATIONS_HPP_INCLUDED
21 #define BOOST_NUMERIC_ODEINT_EXTERNAL_MKL_MKL_OPERATIONS_HPP_INCLUDED
25 #include <mkl_cblas.h>
26 #include <boost/numeric/odeint/algebra/default_operations.hpp>
28 /* exemplary example for writing bindings to the Intel MKL library
29 * see test/mkl for how to use mkl with odeint
30 * this is a quick and dirty implementation showing the general possibility.
31 * It works only with containers based on double and sequential memory allocation.
38 /* only defined for doubles */
41 //template< class Fac1 , class Fac2 > struct scale_sum2;
43 template< class F1 = double , class F2 = F1 >
51 scale_sum2( const Fac1 alpha1 , const Fac2 alpha2 ) : m_alpha1( alpha1 ) , m_alpha2( alpha2 ) { }
53 template< class T1 , class T2 , class T3 >
54 void operator()( T1 &t1 , const T2 &t2 , const T3 &t3) const
55 { // t1 = m_alpha1 * t2 + m_alpha2 * t3;
56 // we get Containers that have size() and [i]-access
58 const int n = t1.size();
59 //boost::numeric::odeint::copy( t1 , t3 );
60 if( &(t2[0]) != &(t1[0]) )
62 cblas_dcopy( n , &(t2[0]) , 1 , &(t1[0]) , 1 );
64 cblas_dscal( n , m_alpha1 , &(t1[0]) , 1 );
65 cblas_daxpy( n , m_alpha2 , &(t3[0]) , 1 , &(t1[0]) , 1 );
66 //daxpby( &n , &m_alpha2 , &(t3[0]) , &one , &m_alpha1 , &(t1[0]) , &one );
70 template< class F1 = double , class F2 = F1 , class F3 = F2 >
80 scale_sum3( const Fac1 alpha1 , const Fac2 alpha2 , const Fac3 alpha3 )
81 : m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) { }
83 template< class T1 , class T2 , class T3 , class T4 >
84 void operator()( T1 &t1 , const T2 &t2 , const T3 &t3 , const T4 &t4 ) const
85 { // t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4;
86 // we get Containers that have size() and [i]-access
88 const int n = t1.size();
89 //boost::numeric::odeint::copy( t1 , t3 );
90 if( &(t2[0]) != &(t1[0]) )
92 cblas_dcopy( n , &(t2[0]) , 1 , &(t1[0]) , 1 );
94 cblas_dscal( n , m_alpha1 , &(t1[0]) , 1 );
95 cblas_daxpy( n , m_alpha2 , &(t3[0]) , 1 , &(t1[0]) , 1 );
96 //daxpby( &n , &m_alpha2 , &(t3[0]) , &one , &m_alpha1 , &(t1[0]) , &one );
97 cblas_daxpy( n , m_alpha3 , &(t4[0]) , 1 , &(t1[0]) , 1 );
101 template< class F1 = double , class F2 = F1 , class F3 = F2 , class F4 = F3 >
113 scale_sum4( const Fac1 alpha1 , const Fac2 alpha2 , const Fac3 alpha3 , const Fac4 alpha4 )
114 : m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) , m_alpha4( alpha4 ) { }
116 template< class T1 , class T2 , class T3 , class T4 , class T5 >
117 void operator()( T1 &t1 , const T2 &t2 , const T3 &t3 , const T4 &t4 , const T5 &t5 ) const
118 { // t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5;
119 // we get Containers that have size() and [i]-access
121 const int n = t1.size();
122 //boost::numeric::odeint::copy( t1 , t3 );
123 if( &(t2[0]) != &(t1[0]) )
125 cblas_dcopy( n , &(t2[0]) , 1 , &(t1[0]) , 1 );
128 cblas_dscal( n , m_alpha1 , &(t1[0]) , 1 );
129 cblas_daxpy( n , m_alpha2 , &(t3[0]) , 1 , &(t1[0]) , 1 );
130 //daxpby( &n , &m_alpha2 , &(t3[0]) , &one , &m_alpha1 , &(t1[0]) , &one );
131 cblas_daxpy( n , m_alpha3 , &(t4[0]) , 1 , &(t1[0]) , 1 );
132 cblas_daxpy( n , m_alpha4 , &(t5[0]) , 1 , &(t1[0]) , 1 );
137 template< class F1 = double , class F2 = F1 , class F3 = F2 , class F4 = F3 , class F5 = F4 >
151 scale_sum5( const Fac1 alpha1 , const Fac2 alpha2 , const Fac3 alpha3 , const Fac4 alpha4 , const Fac5 alpha5 )
152 : m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) , m_alpha4( alpha4 ) , m_alpha5( alpha5 ) { }
154 template< class T1 , class T2 , class T3 , class T4 , class T5 , class T6 >
155 void operator()( T1 &t1 , const T2 &t2 , const T3 &t3 , const T4 &t4 , const T5 &t5 , const T6 &t6 ) const
156 { // t1 = m_alpha1 * t2 + m_alpha2 * t3 + m_alpha3 * t4 + m_alpha4 * t5 + m_alpha5 * t6;
157 // we get Containers that have size() and [i]-access
159 const int n = t1.size();
160 //boost::numeric::odeint::copy( t1 , t3 );
161 if( &(t2[0]) != &(t1[0]) )
163 cblas_dcopy( n , &(t2[0]) , 1 , &(t1[0]) , 1 );
166 cblas_dscal( n , m_alpha1 , &(t1[0]) , 1 );
167 cblas_daxpy( n , m_alpha2 , &(t3[0]) , 1 , &(t1[0]) , 1 );
168 //daxpby( &n , &m_alpha2 , &(t3[0]) , &one , &m_alpha1 , &(t1[0]) , &one );
169 cblas_daxpy( n , m_alpha3 , &(t4[0]) , 1 , &(t1[0]) , 1 );
170 cblas_daxpy( n , m_alpha4 , &(t5[0]) , 1 , &(t1[0]) , 1 );
171 cblas_daxpy( n , m_alpha5 , &(t6[0]) , 1 , &(t1[0]) , 1 );
181 #endif // BOOST_NUMERIC_ODEINT_EXTERNAL_MKL_MKL_OPERATIONS_HPP_INCLUDED