Merge pull request #186 from mmlanger/issue185_bs_performance

Performance improvement for Bulirsch-Stoer-Stepper
headmyshoulder · Dec 28, 2015 · ae53eb1 · ae53eb1
2 parents eb2804b + a625f53
commit ae53eb1
Show file tree

Hide file tree

Showing 2 changed files with 13 additions and 2 deletions.
diff --git a/include/boost/numeric/odeint/stepper/bulirsch_stoer.hpp b/include/boost/numeric/odeint/stepper/bulirsch_stoer.hpp
@@ -98,6 +98,7 @@ class bulirsch_stoer {
           m_interval_sequence( m_k_max+1 ) ,
           m_coeff( m_k_max+1 ) ,
           m_cost( m_k_max+1 ) ,
+          m_facmin_table( m_k_max+1 ) ,
           m_table( m_k_max ) ,
           STEPFAC1( 0.65 ) , STEPFAC2( 0.94 ) , STEPFAC3( 0.02 ) , STEPFAC4( 4.0 ) , KFAC1( 0.8 ) , KFAC2( 0.9 )
     {
@@ -112,6 +113,7 @@ class bulirsch_stoer {
             else
                 m_cost[i] = m_cost[i-1] + m_interval_sequence[i];
             m_coeff[i].resize(i);
+            m_facmin_table[i] = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , 1.0 / static_cast< value_type > ( 2*i+1 ) );
             for( size_t k = 0 ; k < i ; ++k  )
             {
                 const value_type r = static_cast< value_type >( m_interval_sequence[i] ) / static_cast< value_type >( m_interval_sequence[k] );
@@ -415,7 +417,7 @@ class bulirsch_stoer {
         BOOST_USING_STD_MAX();
         using std::pow;
         value_type expo( 1.0/(2*k+1) );
-        value_type facmin = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , expo );
+        value_type facmin = m_facmin_table[k];
         value_type fac;
         if (error == 0.0)
             fac=1.0/facmin;
@@ -504,6 +506,7 @@ class bulirsch_stoer {
     int_vector m_interval_sequence; // stores the successive interval counts
     value_matrix m_coeff;
     int_vector m_cost; // costs for interval count
+    value_vector m_facmin_table; // for precomputed facmin to save pow calls
 
     state_table_type m_table; // sequence of states for extrapolation
 

diff --git a/include/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp b/include/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp
@@ -110,6 +110,7 @@ class bulirsch_stoer_dense_out {
           m_interval_sequence( m_k_max+1 ) ,
           m_coeff( m_k_max+1 ) ,
           m_cost( m_k_max+1 ) ,
+          m_facmin_table( m_k_max+1 ) ,
           m_table( m_k_max ) ,
           m_mp_states( m_k_max+1 ) ,
           m_derivs( m_k_max+1 ) ,
@@ -125,9 +126,15 @@ class bulirsch_stoer_dense_out {
             m_interval_sequence[i] = 2 + 4*i;  // 2 6 10 14 ...
             m_derivs[i].resize( m_interval_sequence[i] );
             if( i == 0 )
+            {
                 m_cost[i] = m_interval_sequence[i];
+                m_facmin_table[0] = 1; // never to be used
+            }
             else
+            {
+                m_facmin_table[i] = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , static_cast<value_type>(1)/(m_interval_sequence[i-1]) );
                 m_cost[i] = m_cost[i-1] + m_interval_sequence[i];
+            }
             m_coeff[i].resize(i);
             for( size_t k = 0 ; k < i ; ++k  )
             {
@@ -429,7 +436,7 @@ class bulirsch_stoer_dense_out {
         using std::pow;
 
         value_type expo = static_cast<value_type>(1)/(m_interval_sequence[k-1]);
-        value_type facmin = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , expo );
+        value_type facmin = m_facmin_table[k];
         value_type fac;
         if (error == 0.0)
             fac = static_cast<value_type>(1)/facmin;
@@ -692,6 +699,7 @@ class bulirsch_stoer_dense_out {
     int_vector m_interval_sequence; // stores the successive interval counts
     value_matrix m_coeff;
     int_vector m_cost; // costs for interval count
+    value_vector m_facmin_table; // for precomputed facmin to save pow calls
 
     state_vector_type m_table; // sequence of states for extrapolation