1 /* 2 [auto_generated] 3 boost/numeric/odeint/external/mkl/mkl_operations.hpp 4 5 [begin_description] 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/ 9 [end_description] 10 11 Copyright 2010-2011 Mario Mulansky 12 Copyright 2011-2013 Karsten Ahnert 13 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) 17 */ 18 19 20 #ifndef BOOST_NUMERIC_ODEINT_EXTERNAL_MKL_MKL_OPERATIONS_HPP_INCLUDED 21 #define BOOST_NUMERIC_ODEINT_EXTERNAL_MKL_MKL_OPERATIONS_HPP_INCLUDED 22 23 #include <iostream> 24 25 #include <mkl_cblas.h> 26 #include <boost/numeric/odeint/algebra/default_operations.hpp> 27 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. 32 */ 33 34 namespace boost { 35 namespace numeric { 36 namespace odeint { 37 38 /* only defined for doubles */ 39 struct mkl_operations 40 { 41 //template< class Fac1 , class Fac2 > struct scale_sum2; 42 43 template< class F1 = double , class F2 = F1 > 44 struct scale_sum2 45 { 46 typedef double Fac1; 47 typedef double Fac2; 48 const Fac1 m_alpha1; 49 const Fac2 m_alpha2; 50 scale_sum2boost::numeric::odeint::mkl_operations::scale_sum251 scale_sum2( const Fac1 alpha1 , const Fac2 alpha2 ) : m_alpha1( alpha1 ) , m_alpha2( alpha2 ) { } 52 53 template< class T1 , class T2 , class T3 > operator ()boost::numeric::odeint::mkl_operations::scale_sum254 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 57 58 const int n = t1.size(); 59 //boost::numeric::odeint::copy( t1 , t3 ); 60 if( &(t2[0]) != &(t1[0]) ) 61 { 62 cblas_dcopy( n , &(t2[0]) , 1 , &(t1[0]) , 1 ); 63 } 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 ); 67 } 68 }; 69 70 template< class F1 = double , class F2 = F1 , class F3 = F2 > 71 struct scale_sum3 72 { 73 typedef double Fac1; 74 typedef double Fac2; 75 typedef double Fac3; 76 const Fac1 m_alpha1; 77 const Fac2 m_alpha2; 78 const Fac3 m_alpha3; 79 scale_sum3boost::numeric::odeint::mkl_operations::scale_sum380 scale_sum3( const Fac1 alpha1 , const Fac2 alpha2 , const Fac3 alpha3 ) 81 : m_alpha1( alpha1 ) , m_alpha2( alpha2 ) , m_alpha3( alpha3 ) { } 82 83 template< class T1 , class T2 , class T3 , class T4 > operator ()boost::numeric::odeint::mkl_operations::scale_sum384 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 87 88 const int n = t1.size(); 89 //boost::numeric::odeint::copy( t1 , t3 ); 90 if( &(t2[0]) != &(t1[0]) ) 91 { 92 cblas_dcopy( n , &(t2[0]) , 1 , &(t1[0]) , 1 ); 93 } 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 ); 98 } 99 }; 100 101 template< class F1 = double , class F2 = F1 , class F3 = F2 , class F4 = F3 > 102 struct scale_sum4 103 { 104 typedef double Fac1; 105 typedef double Fac2; 106 typedef double Fac3; 107 typedef double Fac4; 108 const Fac1 m_alpha1; 109 const Fac2 m_alpha2; 110 const Fac3 m_alpha3; 111 const Fac4 m_alpha4; 112 scale_sum4boost::numeric::odeint::mkl_operations::scale_sum4113 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 ) { } 115 116 template< class T1 , class T2 , class T3 , class T4 , class T5 > operator ()boost::numeric::odeint::mkl_operations::scale_sum4117 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 120 121 const int n = t1.size(); 122 //boost::numeric::odeint::copy( t1 , t3 ); 123 if( &(t2[0]) != &(t1[0]) ) 124 { 125 cblas_dcopy( n , &(t2[0]) , 1 , &(t1[0]) , 1 ); 126 } 127 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 ); 133 } 134 }; 135 136 137 template< class F1 = double , class F2 = F1 , class F3 = F2 , class F4 = F3 , class F5 = F4 > 138 struct scale_sum5 139 { 140 typedef double Fac1; 141 typedef double Fac2; 142 typedef double Fac3; 143 typedef double Fac4; 144 typedef double Fac5; 145 const Fac1 m_alpha1; 146 const Fac2 m_alpha2; 147 const Fac3 m_alpha3; 148 const Fac4 m_alpha4; 149 const Fac5 m_alpha5; 150 scale_sum5boost::numeric::odeint::mkl_operations::scale_sum5151 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 ) { } 153 154 template< class T1 , class T2 , class T3 , class T4 , class T5 , class T6 > operator ()boost::numeric::odeint::mkl_operations::scale_sum5155 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 158 159 const int n = t1.size(); 160 //boost::numeric::odeint::copy( t1 , t3 ); 161 if( &(t2[0]) != &(t1[0]) ) 162 { 163 cblas_dcopy( n , &(t2[0]) , 1 , &(t1[0]) , 1 ); 164 } 165 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 ); 172 } 173 }; 174 175 }; 176 177 } // odeint 178 } // numeric 179 } // boost 180 181 #endif // BOOST_NUMERIC_ODEINT_EXTERNAL_MKL_MKL_OPERATIONS_HPP_INCLUDED 182