runge_kutta.h

Go to the documentation of this file.

#ifndef _DG_RK_
#define _DG_RK_
 
#include <cassert>
#include <array>
#include <tuple>
 
#include "ode.h"
#include "backend/exceptions.h"
#include "tableau.h"
#include "blas1.h"
#include "implicit.h"
 
namespace dg{
namespace detail{
// this is a dense matrix-matrix multiplication
// a generalization of the corresponding blas2::gemv DenseMatrix algorithm
template<class ContainerType>
void gemm(
        const std::array<double,2>& alpha,
        const DenseMatrix<ContainerType>& m,
        const std::array<const std::vector<double>*,2> x,
        const std::array<double,2>& beta,
        std::array<ContainerType*,2>&& y)
{
    const unsigned size = m.num_cols();
    unsigned i=0;
    if( size >= 8)
        for( i=0; i<size/8; i++)
            dg::blas1::subroutine( dg::EmbeddedPairSum(),
                    *y[0], *y[1], i==0 ? beta[0] : 1., i==0 ? beta[1] : 1.,
                    alpha[0]*(*x[0])[i*8+0], alpha[1]*(*x[1])[i*8+0], *m.get()[i*8+0],
                    alpha[0]*(*x[0])[i*8+1], alpha[1]*(*x[1])[i*8+1], *m.get()[i*8+1],
                    alpha[0]*(*x[0])[i*8+2], alpha[1]*(*x[1])[i*8+2], *m.get()[i*8+2],
                    alpha[0]*(*x[0])[i*8+3], alpha[1]*(*x[1])[i*8+3], *m.get()[i*8+3],
                    alpha[0]*(*x[0])[i*8+4], alpha[1]*(*x[1])[i*8+4], *m.get()[i*8+4],
                    alpha[0]*(*x[0])[i*8+5], alpha[1]*(*x[1])[i*8+5], *m.get()[i*8+5],
                    alpha[0]*(*x[0])[i*8+6], alpha[1]*(*x[1])[i*8+6], *m.get()[i*8+6],
                    alpha[0]*(*x[0])[i*8+7], alpha[1]*(*x[1])[i*8+7], *m.get()[i*8+7]);
    unsigned l=0;
    if( size%8 >= 4)
        for( l=0; l<(size%8)/4; l++)
            dg::blas1::subroutine( dg::EmbeddedPairSum(),
                    *y[0], *y[1], size < 8 ?  beta[0]  : 1., size < 8 ? beta[1] : 1.,
                    alpha[0]*(*x[0])[i*8+l*4+0], alpha[1]*(*x[1])[i*8+l*4+0], *m.get()[i*8+l*4+0],
                    alpha[0]*(*x[0])[i*8+l*4+1], alpha[1]*(*x[1])[i*8+l*4+1], *m.get()[i*8+l*4+1],
                    alpha[0]*(*x[0])[i*8+l*4+2], alpha[1]*(*x[1])[i*8+l*4+2], *m.get()[i*8+l*4+2],
                    alpha[0]*(*x[0])[i*8+l*4+3], alpha[1]*(*x[1])[i*8+l*4+3], *m.get()[i*8+l*4+3]);
    unsigned k=0;
    if( (size%8)%4 >= 2)
        for( k=0; k<((size%8)%4)/2; k++)
            dg::blas1::subroutine( dg::EmbeddedPairSum(),
                    *y[0], *y[1], size < 4 ?  beta[0]  : 1., size < 4 ? beta[1] : 1.,
                    alpha[0]*(*x[0])[i*8+l*4+k*2+0], alpha[1]*(*x[1])[i*8+l*4+k*2+0], *m.get()[i*8+l*4+k*2+0],
                    alpha[0]*(*x[0])[i*8+l*4+k*2+1], alpha[1]*(*x[1])[i*8+l*4+k*2+1], *m.get()[i*8+l*4+k*2+1]);
    if( ((size%8)%4)%2 == 1)
        dg::blas1::subroutine( dg::EmbeddedPairSum(),
                    *y[0], *y[1], size < 2 ?  beta[0]  : 1., size < 2 ? beta[1] : 1.,
                    alpha[0]*(*x[0])[i*8+l*4+k*2], alpha[1]*(*x[1])[i*8+l*4+k*2], *m.get()[i*8+l*4+k*2]);
 
}
} // namespace detail
 
struct IdentityFilter
{
    template<class ContainerType1>
    void operator()( ContainerType1&) const{ }
};
template<class ContainerType>
struct FilteredERKStep;
//Should we try if filters inside RHS are equally usable?
 
 
 
template< class ContainerType>
struct ERKStep
{
    using value_type = get_value_type<ContainerType>;
    using container_type = ContainerType; 
    ERKStep() = default;
    ERKStep( ConvertsToButcherTableau<value_type> tableau, const ContainerType&
        copyable): m_ferk(tableau, copyable)
    {
    }
 
    template<class ...Params>
    void construct( Params&& ...ps)
    {
        //construct and swap
        *this = ERKStep( std::forward<Params>( ps)...);
    }
    const ContainerType& copyable()const{ return m_ferk.copyable();}
 
    void ignore_fsal(){ m_ferk.ignore_fsal();}
    void enable_fsal(){ m_ferk.enable_fsal();}
 
    template<class ExplicitRHS>
    void step( ExplicitRHS& rhs, value_type t0, const ContainerType& u0, value_type& t1, ContainerType& u1, value_type dt, ContainerType& delta)
    {
        dg::IdentityFilter id;
        m_ferk.step( std::tie( rhs, id), t0, u0, t1, u1, dt, delta);
    }
    template<class ExplicitRHS>
    void step( ExplicitRHS& rhs, value_type t0, const ContainerType& u0, value_type&
            t1, ContainerType& u1, value_type dt)
    {
        dg::IdentityFilter id;
        m_ferk.step( std::tie( rhs, id), t0, u0, t1, u1, dt);
    }
    unsigned order() const {
        return m_ferk.order();
    }
    unsigned embedded_order() const {
        return m_ferk.embedded_order();
    }
    unsigned num_stages() const{
        return m_ferk.num_stages();
    }
  private:
    FilteredERKStep<ContainerType> m_ferk;
};
 
template< class ContainerType>
struct FilteredERKStep
{
    using value_type = get_value_type<ContainerType>;
    using container_type = ContainerType; 
    FilteredERKStep() { m_k.resize(1); }
    FilteredERKStep( ConvertsToButcherTableau<value_type> tableau, const
            ContainerType& copyable): m_rk(tableau), m_k(m_rk.num_stages(),
                copyable)
    {
        m_rkb.resize(m_k.size()), m_rkd.resize( m_k.size());
        for( unsigned i=0; i<m_k.size(); i++)
        {
            m_rkb[i] = m_rk.b(i);
            m_rkd[i] = m_rk.d(i);
        }
    }
 
    template<class ...Params>
    void construct( Params&& ...ps)
    {
        //construct and swap
        *this = FilteredERKStep( std::forward<Params>( ps)...);
    }
    const ContainerType& copyable()const{ return m_k[0];}
 
    void ignore_fsal(){ m_ignore_fsal = true;}
    void enable_fsal(){ m_ignore_fsal = false;}
 
    template<class ExplicitRHS, class Limiter>
    void step( const std::tuple<ExplicitRHS,Limiter>& rhs, value_type t0, const ContainerType& u0, value_type& t1, ContainerType& u1, value_type dt, ContainerType& delta)
    {
        step ( rhs, t0, u0, t1, u1, dt, delta, true);
    }
    template<class ExplicitRHS, class Limiter>
    void step( const std::tuple<ExplicitRHS, Limiter>& rhs, value_type t0, const ContainerType& u0, value_type&
            t1, ContainerType& u1, value_type dt)
    {
        if( !m_tmp_allocated)
        {
            dg::assign( m_k[0], m_tmp);
            m_tmp_allocated = true;
        }
        step ( rhs, t0, u0, t1, u1, dt, m_tmp, false);
    }
    unsigned order() const {
        return m_rk.order();
    }
    unsigned embedded_order() const {
        return m_rk.embedded_order();
    }
    unsigned num_stages() const{
        return m_rk.num_stages();
    }
  private:
    template<class ExplicitRHS, class Limiter>
    void step( const std::tuple<ExplicitRHS, Limiter>& rhs, value_type t0, const ContainerType& u0, value_type& t1, ContainerType& u1, value_type dt, ContainerType& delta, bool);
    ButcherTableau<value_type> m_rk;
    std::vector<value_type> m_rkb, m_rkd;
    std::vector<ContainerType> m_k;
    value_type m_t1 = 1e300;//remember the last timestep at which ERK is called
    bool m_ignore_fsal = false;
    ContainerType m_tmp; //  only conditionally allocated
    bool m_tmp_allocated = false;
};
 
 
template< class ContainerType>
template<class ExplicitRHS, class Limiter>
void FilteredERKStep<ContainerType>::step( const std::tuple<ExplicitRHS, Limiter>& ode, value_type t0, const ContainerType& u0, value_type& t1, ContainerType& u1, value_type dt, ContainerType& delta, bool compute_delta)
{
    unsigned s = m_rk.num_stages();
    std::vector<const ContainerType*> k_ptrs = dg::asPointers( m_k);
    //0 stage: probe
    value_type tu = t0;
    if( t0 != m_t1 || m_ignore_fsal)
        std::get<0>(ode)(t0, u0, m_k[0]); //freshly compute k_0
    //else take from last call
    for ( unsigned i=1; i<s; i++)
    {
        std::vector<value_type> rka( i);
        for( unsigned l=0; l<i; l++)
            rka[l] = m_rk.a(i,l);
 
        tu = DG_FMA( dt,m_rk.c(i),t0); //l=0
        dg::blas1::copy( u0, delta); // can't use u1 here cause u0 can alias
        dg::blas2::gemv( dt, dg::asDenseMatrix( k_ptrs, i), rka, 1., delta);
 
        std::get<1>(ode)( delta);
        std::get<0>(ode)( tu, delta, m_k[i]);
    }
    //Now add everything up to get solution and error estimate
    dg::blas1::copy( u0, u1);
    if( compute_delta)
        detail::gemm( {dt,dt}, dg::asDenseMatrix(k_ptrs), {&m_rkb, &m_rkd},
            {1.,0.}, {&u1, &delta});
    else
        blas2::gemv( dt, dg::asDenseMatrix(k_ptrs), m_rkb, 1., u1);
    std::get<1>(ode)( u1);
    //make sure (t1,u1) is the last call to f
    m_t1 = t1 = t0 + dt;
    if(!m_rk.isFsal() )
        std::get<0>(ode)( t1, u1, m_k[0]);
    else
    {
        using std::swap;
        swap( m_k[0], m_k[s-1]); //enable free swap functions
    }
}
 
 
template<class ContainerType>
struct ARKStep
{
    using value_type = get_value_type<ContainerType>;
    using container_type = ContainerType; 
    ARKStep(){ m_kI.resize(1); }
    ARKStep( std::string name, const ContainerType& copyable)
    {
        std::string exp_name = name+" (explicit)";
        std::string imp_name = name+" (implicit)";
        *this = ARKStep( exp_name, imp_name, copyable);
    }
    ARKStep( ConvertsToButcherTableau<value_type> ex_tableau,
             ConvertsToButcherTableau<value_type> im_tableau,
             const ContainerType& copyable
             ):
         m_rkE(ex_tableau),
         m_rkI(im_tableau),
         m_kE(m_rkE.num_stages(), copyable),
         m_kI(m_rkI.num_stages(), copyable)
    {
        assert( m_rkE.num_stages() == m_rkI.num_stages());
        // check fsal
        assert( m_rkI.a(0,0) == 0);
        assert( m_rkI.c(m_rkI.num_stages()-1) == 1);
        check_implicit_fsal();
        assign_coeffs();
    }
    template<class ...Params>
    void construct( Params&& ...ps)
    {
        //construct and swap
        *this = ARKStep( std::forward<Params>( ps)...);
    }
    const ContainerType& copyable()const{ return m_kI[0];}
 
    template< class ExplicitRHS, class ImplicitRHS, class Solver>
    void step( const std::tuple<ExplicitRHS, ImplicitRHS, Solver>& ode, value_type t0, const ContainerType& u0, value_type& t1, ContainerType& u1, value_type dt, ContainerType& delta);
    unsigned order() const {
        return m_rkE.order();
    }
    unsigned embedded_order() const {
        return m_rkE.order();
    }
    unsigned num_stages() const{
        return m_rkE.num_stages();
    }
    private:
    ButcherTableau<value_type> m_rkE, m_rkI;
    std::vector<ContainerType> m_kE, m_kI;
    std::vector<value_type> m_rkb, m_rkd;
    value_type m_t1 = 1e300;
    bool m_implicit_fsal = false;
    void check_implicit_fsal(){
        m_implicit_fsal = true;
        for( unsigned i=0; i<m_rkI.num_stages(); i++)
            if( m_rkI.a(i,0) != 0)
                m_implicit_fsal = false;
    }
    void assign_coeffs()
    {
        m_rkb.resize( 2*m_rkI.num_stages());
        m_rkd.resize( 2*m_rkI.num_stages());
        for( unsigned i=0; i<m_rkI.num_stages(); i++)
        {
            m_rkb[2*i] = m_rkE.b(i);
            m_rkb[2*i+1] = m_rkI.b(i);
            m_rkd[2*i] = m_rkE.d(i);
            m_rkd[2*i+1] = m_rkI.d(i);
        }
    }
};
 
template<class ContainerType>
template< class Explicit, class Implicit, class Solver>
void ARKStep<ContainerType>::step( const std::tuple<Explicit,Implicit,Solver>& ode, value_type t0, const ContainerType& u0, value_type& t1, ContainerType& u1, value_type dt, ContainerType& delta)
{
    unsigned s = m_rkE.num_stages();
    value_type tu = t0;
    //0 stage
    if( t0 != m_t1)
        std::get<0>(ode)(t0, u0, m_kE[0]); //freshly compute k_0
    if( !m_implicit_fsal) // all a(i,0) == 0
        std::get<1>(ode)(t0, u0, m_kI[0]);
    // DO NOT HOLD POINTERS AS PRVIATE MEMBERS
    std::vector<const ContainerType*> k_ptrs( 2*m_rkI.num_stages());
    for( unsigned i=0; i<m_rkI.num_stages(); i++)
    {
        k_ptrs[2*i  ] = &m_kE[i];
        k_ptrs[2*i+1] = &m_kI[i];
    }
 
    for( unsigned i=1; i<s; i++)
    {
        std::vector<value_type> rka( 2*i);
        for(unsigned l=0; l<i; l++)
        {
            rka[2*l]   = m_rkE.a(i,l);
            rka[2*l+1] = m_rkI.a(i,l);
        }
        tu = DG_FMA( m_rkI.c(i),dt, t0);
        dg::blas1::copy( u0, m_kI[i]);
        dg::blas2::gemv( dt, dg::asDenseMatrix( k_ptrs, 2*i), rka, 1., m_kI[i]);
        value_type alpha = dt*m_rkI.a(i,i);
        std::get<2>(ode)( alpha, tu, delta, m_kI[i]);
        dg::blas1::axpby( 1./alpha, delta, -1./alpha, m_kI[i]);
        std::get<0>(ode)(tu, delta, m_kE[i]);
    }
    m_t1 = t1 = tu;
    //Now compute result and error estimate
 
    dg::blas1::copy( u0, u1);
    detail::gemm( {dt,dt}, dg::asDenseMatrix(k_ptrs), {&m_rkb, &m_rkd},
            {1.,0.}, {&u1, &delta});
    //make sure (t1,u1) is the last call to ex
    std::get<0>(ode)(t1,u1,m_kE[0]);
}
 
template<class ContainerType>
struct DIRKStep
{
    //MW Dirk methods cannot have stage order greater than 1
    using value_type = get_value_type<ContainerType>;
    using container_type = ContainerType; 
    DIRKStep(){ m_kI.resize(1); }
 
    DIRKStep( ConvertsToButcherTableau<value_type> im_tableau,
               const ContainerType& copyable):
         m_rkI(im_tableau),
         m_kI(m_rkI.num_stages(), copyable)
    {
        m_rkIb.resize(m_kI.size()), m_rkId.resize(m_kI.size());
        for( unsigned i=0; i<m_kI.size(); i++)
        {
            m_rkIb[i] = m_rkI.b(i);
            m_rkId[i] = m_rkI.d(i);
        }
    }
 
    template<class ...Params>
    void construct( Params&& ...ps)
    {
        //construct and swap
        *this = DIRKStep( std::forward<Params>( ps)...);
    }
    const ContainerType& copyable()const{ return m_kI[0];}
 
    template< class ImplicitRHS, class Solver>
    void step( const std::tuple<ImplicitRHS,Solver>& ode, value_type t0, const
        ContainerType& u0, value_type& t1, ContainerType& u1, value_type dt,
        ContainerType& delta)
    {
        step( ode, t0, u0, t1, u1, dt, delta, true);
    }
    template< class ImplicitRHS, class Solver>
    void step( const std::tuple<ImplicitRHS, Solver>& ode, value_type t0, const
        ContainerType& u0, value_type& t1, ContainerType& u1, value_type dt)
    {
        if( !m_allocated)
        {
            dg::assign( m_kI[0], m_tmp);
            m_allocated = true;
        }
        step( ode, t0, u0, t1, u1, dt, m_tmp, false);
    }
    unsigned order() const {
        return m_rkI.order();
    }
    unsigned embedded_order() const {
        return m_rkI.order();
    }
    unsigned num_stages() const{
        return m_rkI.num_stages();
    }
 
    private:
    template< class ImplicitRHS, class Solver>
    void step( const std::tuple<ImplicitRHS, Solver>& ode, value_type t0, const
            ContainerType& u0, value_type& t1, ContainerType& u1, value_type
            dt, ContainerType& delta, bool compute_delta);
    ButcherTableau<value_type> m_rkI;
    std::vector<ContainerType> m_kI;
    ContainerType m_tmp;
    std::vector<value_type> m_rkIb, m_rkId;
    bool m_allocated = false;
};
 
template<class ContainerType>
template< class ImplicitRHS, class Solver>
void DIRKStep<ContainerType>::step( const std::tuple<ImplicitRHS,Solver>& ode,  value_type t0, const ContainerType& u0, value_type& t1, ContainerType& u1, value_type dt, ContainerType& delta, bool compute_delta)
{
    unsigned s = m_rkI.num_stages();
    value_type tu = t0;
    //0 stage
    //rhs = u0
    tu = DG_FMA( m_rkI.c(0),dt, t0);
    value_type alpha;
    if( m_rkI.a(0,0) !=0 )
    {
        alpha = dt*m_rkI.a(0,0);
        std::get<1>(ode)( alpha, tu, delta, u0);
        dg::blas1::axpby( 1./alpha, delta, -1./alpha, u0, m_kI[0]);
    }
    else
        std::get<0>(ode)(tu, u0, m_kI[0]);
    std::vector<const ContainerType*> kIptr = dg::asPointers( m_kI);
 
    for( unsigned i=1; i<s; i++)
    {
        tu = DG_FMA( m_rkI.c(i),dt, t0);
        dg::blas1::copy( u0, m_kI[i]);
        std::vector<value_type> rkIa( i);
        for( unsigned l=0; l<i; l++)
            rkIa[l] = m_rkI.a(i,l);
        dg::blas2::gemv( dt, dg::asDenseMatrix(kIptr,i), rkIa, 1., m_kI[i]);
        if( m_rkI.a(i,i) !=0 )
        {
            alpha = dt*m_rkI.a(i,i);
            std::get<1>(ode)( alpha, tu, delta, m_kI[i]);
            dg::blas1::axpby( 1./alpha, delta, -1./alpha, m_kI[i]);
        }
        else
            std::get<0>(ode)(tu, delta, m_kI[i]);
    }
    t1 = t0 + dt;
    //Now compute result and error estimate
    dg::blas1::copy( u0, u1);
    if( compute_delta)
        detail::gemm( {dt,dt}, dg::asDenseMatrix(kIptr), {&m_rkIb, &m_rkId},
            {1.,0.}, {&u1, &delta});
    else
        blas2::gemv( dt, dg::asDenseMatrix(kIptr), m_rkIb, 1., u1);
}
 
template<class ContainerType>
using RungeKutta = ERKStep<ContainerType>;
 
template<class ContainerType>
using FilteredRungeKutta = FilteredERKStep<ContainerType>;
 
template<class ContainerType>
struct ShuOsher
{
    using value_type = get_value_type<ContainerType>;
    using container_type = ContainerType; 
    ShuOsher(){m_u.resize(1);}
    ShuOsher( dg::ConvertsToShuOsherTableau<value_type> tableau, const ContainerType& copyable): m_t( tableau), m_u(  m_t.num_stages(), copyable), m_k(m_u)
        { }
    template<class ...Params>
    void construct( Params&& ...ps)
    {
        //construct and swap
        *this = ShuOsher( std::forward<Params>( ps)...);
    }
    const ContainerType& copyable()const{ return m_u[0];}
 
    template<class ExplicitRHS, class Limiter>
    void step( const std::tuple<ExplicitRHS, Limiter>& ode, value_type t0, const ContainerType& u0, value_type& t1, ContainerType& u1, value_type dt){
        unsigned s = m_t.num_stages();
        std::vector<value_type> ts( m_t.num_stages()+1);
        ts[0] = t0;
        dg::blas1::copy(u0, m_u[0]);
        if( t0 != m_t1 ) //this is the first time we call step
            std::get<0>(ode)(ts[0], m_u[0], m_k[0]); //freshly compute k_0
        for( unsigned i=1; i<=s; i++)
        {
            dg::blas1::axpbypgz( m_t.alpha(i-1,0), m_u[0], dt*m_t.beta(i-1,0), m_k[0], 0., i==s ? u1 : m_u[i]);
            ts[i] = m_t.alpha(i-1,0)*ts[0] + dt*m_t.beta(i-1,0);
            for( unsigned j=1; j<i; j++)
            {
                //about the i-1: it is unclear to me how the ShuOsher tableau makes implicit schemes
                dg::blas1::axpbypgz( m_t.alpha(i-1,j), m_u[j], dt*m_t.beta(i-1,j), m_k[j], 1., i==s ? u1 : m_u[i]);
                ts[i] += m_t.alpha(i-1,j)*ts[j] + dt*m_t.beta(i-1,j);
 
            }
            std::get<1>(ode)( i==s ? u1 : m_u[i]);
            if(i!=s)
                std::get<0>(ode)(ts[i], m_u[i], m_k[i]);
            else
                //make sure (t1,u1) is the last call to f
                std::get<0>(ode)(ts[i], u1, m_k[0]);
        }
        m_t1 = t1 = ts[s];
    }
    unsigned order() const {
        return m_t.order();
    }
    unsigned num_stages() const{
        return m_t.num_stages();
    }
  private:
    ShuOsherTableau<value_type> m_t;
    std::vector<ContainerType> m_u, m_k;
    value_type m_t1 = 1e300;
};
template<class ContainerType>
using ImplicitRungeKutta = DIRKStep<ContainerType>;
 
 
inline bool is_same( double x, double y, double eps = 1e-15)
{
    return fabs(x - y) < eps * std::max(1.0, std::max( fabs(x), fabs(y)));
}
inline bool is_same( float x, float y, float eps = 1e-6)
{
    return fabsf(x - y) < eps * std::max(1.0f, std::max( fabsf(x), fabsf(y)));
}
 
template<class T>
inline bool is_same( T x, T y)
{
    return x == y;
}
 
inline bool is_divisable( double a, double b, double eps = 1e-15)
{
    return is_same( round(a/b)*b, a, eps);
}
inline bool is_divisable( float a, float b, float eps = 1e-6)
{
    return is_same( (float)round(a/b)*b, (float)a, eps);
}
 
template<class ContainerType>
struct SinglestepTimeloop : public aTimeloop<ContainerType>
{
    using container_type = ContainerType;
    using value_type = dg::get_value_type<ContainerType>;
    SinglestepTimeloop( ) = default;
 
    SinglestepTimeloop( std::function<void ( value_type, const ContainerType&,
                value_type&, ContainerType&, value_type)> step, value_type dt = 0 )
        : m_step(step), m_dt(dt){}
 
    template<class Stepper, class ODE>
    SinglestepTimeloop( Stepper&& stepper, ODE&& ode, value_type dt = 0)
    {
        // Open/Close Principle (OCP) Software entities should be open for extension but closed for modification
        m_step = [=, cap = std::tuple<Stepper, ODE>(std::forward<Stepper>(stepper),
                std::forward<ODE>(ode))  ]( auto t0, const auto& y0, auto& t1, auto& y1,
                auto dtt) mutable
        {
            std::get<0>(cap).step( std::get<1>(cap), t0, y0, t1, y1, dtt);
        };
        m_dt = dt;
    }
 
    template<class ...Params>
    void construct( Params&& ...ps)
    {
        //construct and swap
        *this = SinglestepTimeloop( std::forward<Params>( ps)...);
    }
 
    void set_dt( value_type dt){ m_dt = dt;}
 
    void integrate_steps( value_type t0, const container_type& u0, value_type t1,
            container_type& u1, unsigned steps)
    {
        set_dt( (t1-t0)/(value_type)steps );
        this->integrate( t0, u0, t1, u1);
    }
 
    virtual SinglestepTimeloop* clone() const{ return new
        SinglestepTimeloop(*this);}
    private:
    virtual void do_integrate(value_type& t0, const container_type& u0, value_type t1, container_type& u1, enum to mode) const;
    std::function<void ( value_type, const ContainerType&, value_type&,
            ContainerType&, value_type)> m_step;
    virtual value_type do_dt( ) const { return m_dt;}
    value_type m_dt;
};
 
template< class ContainerType>
void SinglestepTimeloop<ContainerType>::do_integrate(
        value_type&  t_begin, const ContainerType&
        begin, value_type t_end, ContainerType& end,
        enum to mode ) const
{
    bool forward = (t_end - t_begin > 0);
    if( (m_dt < 0 && forward) || ( m_dt > 0 && !forward) )
        throw dg::Error( dg::Message(_ping_)<<"Timestep has wrong sign! dt "<<m_dt);
    dg::blas1::copy( begin, end);
    if( m_dt == 0)
        throw dg::Error( dg::Message(_ping_)<<"Timestep may not be zero in SinglestepTimeloop!");
    if( is_divisable( t_end-t_begin, m_dt))
    {
        unsigned N = (unsigned)round((t_end - t_begin)/m_dt);
        for( unsigned i=0; i<N; i++)
            m_step( t_begin, end, t_begin, end, m_dt);
    }
    else
    {
        unsigned N = (unsigned)floor( (t_end-t_begin)/m_dt);
        for( unsigned i=0; i<N; i++)
            m_step( t_begin, end, t_begin, end, m_dt);
        if( dg::to::exact == mode)
        {
            value_type dt_final = t_end - t_begin;
            m_step( t_begin, end, t_begin, end, dt_final);
        }
        else
            m_step( t_begin, end, t_begin, end, m_dt);
    }
    return;
}
 
 
} //namespace dg
 
#endif //_DG_RK_