Segmentation fault in Bulirsch-Stoer dense output stepper

[ds283]

I am seeing an intermittent segmentation fault in the Bulirsch-Stoer dense output stepper on Linux and OS X.

The segmentation fault occurs in bulirsch_stoer_dense_out::calculate_finite_difference. I believe the intermittency is because the segmentation fault occurs only when the stepper switches to the highest available order m_k_max.

The culprit seems to be operations like this

            m_algebra.for_each2( m_diffs[kappa][j_diffs].m_v , m_derivs[j][m+kappa].m_v ,
                                 typename operations_type::template scale_sum1< value_type >( fac ) );

when kappa takes its maximum value kappa=2*k+1=17 for k=8. I can check that m_diffs.size()=17 and therefore the subscript is out of bounds.

In the constructor for bulirsch_stoer_dense_out, m_diffs is allocated only 2*m_k_max+1 components whereas I believe it should be 2*m_k_max+2. Likewise, the allocation loop at the end of the constructor should presumably be adjusted to read (ie. the start value moved to 2*(m_k_max)+1)

        for( int i = 2*(m_k_max)+1 ; i >=0  ; i-- )
        {
            m_diffs[i].resize( num );
            num += (i+1)%2;
        }

and I assume also in bulirsch_stoer_dense_out:::resize_impl the last for-loop should be (ie. the upper limit of the outer loop moved to 2*m_k_max+2)

        for( size_t i = 0 ; i < 2*m_k_max+2 ; ++i )
            for( size_t j = 0 ; j < m_diffs[i].size() ; ++j )
                resized |= adjust_size_by_resizeability( m_diffs[i][j] , x , typename is_resizeable<deriv_type>::type() );

These changes work for me so far. If they appear to be correct I can submit them as a patch.

[mariomulansky]

Yes. A short look into the code also convinced me that m_diffs should be initialized with 2*(m_k_max+1) instead of 2*m_k_max+1. The other changes then follow immediately. Do you have some minimal example that we could somehow include as test for this bug? Thanks a lot for finding it!

[ds283]

Unfortunately, no - I don't have a minimal example of this kind. The system I'm solving is large and certain parameters trigger the algorithm to go to order m_k_max=8, but I haven't managed to deduce exactly what causes this behaviour.

If the out-of-bounds problem hasn't often been encountered I'm assuming that ODEs which send the algorithm to order m_k_max=8 are not common; but if I sort out the cause I will try to produce a small system which reproduces it.

I've generated a pull request with the fixes.

[mariomulansky]

My guess would be that the switch to maximal order could be triggered by non-continuous points or maybe jumps in the rhs equations. Does your model might lead to such situation for some parameter values?

[ds283]

There aren't explicit jumps, but there is a region where the RHS depends on very accurate cancellations. If the cancellation is imperfect then it could resemble discontinuities.

If I have time I may try to cook up something along these lines.

[mariomulansky]

Great. thanks a lot for your efforts.

[mariomulansky]

Fixed since 656e146