Numerik und Template-Metaprogrammierung

[headmyshoulder]

Einführung

In diesem Artikel werde ich einige Techniken im Zusammenhang mit Templatemetaprogrammierung und numerischen Verfahren vorstellen. Insbesondere wird es dabei um numerische Verfahren zum Lösen von Differentialgleichungen (DGL) gehen. Die hier vorgestellten Methoden lassen sich aber problemlos auf eine Vielzahl von anderen numerischen Problemen anwenden.

Templatemetaprogrammierung ist eine moderne Programmiertechnik, bei welcher der C++ Compiler elementare Berechnungen durchführt oder komplexen C++-Code und Algorithmen erzeugt. Ein entscheidender Aspekt dabei ist, dass Templatemetaprogramme während des Kompiliervorganges ausgewertet bzw. ausgeführt werden. Das klassische Hello-World-Programm ist hier die Berechnung der Fakultät, welche direkt vom Compiler berechnet wird. Zur Laufzeit des fertigen Programmes wird diese Berechnung keine Zeit in Anspruch nehmen.

Templatemetaprogrammierung kann auch für komplexe Probleme und Algorithmen verwendet werden. Ein Beispiel ist die Bibliothek Boost.MetaStateMachine, mit welcher Zustandsmachinen anhand einer (statischen) Übergangstabelle während des Kompilierens erzeugt werden können. Ein weiteres Beispiel ist Boost.Accumulators, eine Bibliothek für statistische Auswertungen. Dort werden die genauen Berechnung und statistischen Routinen auch während der Kompilierphase erzeugt. Außerdem werden Techniken der Templatemetaprogrammierung auch oft im Zusammenhang mit Typkonvertierungen und Typüberprüfungen benutzt.

Es gibt eine Vielzahl von Bibliotheken, welche mehr oder weniger mit dem Thema Templatemetaprogrammierung zu tun haben. Die meisten davon sind in den Boost-Bibliotheken zu finden. Bevor man anfängt Metaprogramme selbst zu schreiben und zu implementieren, empfehle ich aber zuerst dort nachzusehen, ob nicht bereits eine adäquate Lösung existiert. Selbst durch das konsequente Anwenden und Benutzen dieser Bibliotheken kann man sehr viel lernen. Als wichtigste Bibliothek ist dabei sicherlich Boost.MPL mit einer Vielzahl von Metafunktionen, -algorithmen und Compile-Time-Sequenzen zu nennen. Eine weitere wichtige Bibliothek ist Boost.Fusion, in welcher die Runtime-Aspekte mehr betont werden. Zusätzlich müssen in der Templatemetaprogrammierung sehr oft Typen klassifiziert werden. Dafür stehen etliche Funktion in Boost.TypeTraits zur Verfügung. Beispielsweise werden Metafunktionen wie is_pointer , is_class , aber auch Konversionen wie add_const und remove_const zur Verfügung gestellt. Ein Teil dieser Funktionalität wird auch Einzug in den kommenden C++-Standard finden.

Zur Entwicklung von Meta-Algorithmen werden wir hier ausschliesslich Boost.MPL verwenden. Diese Bibliothek ist in ihrem Umfang und ihren Anwendungsmöglichkeiten am weitesten entwickelt. Es gibt schlicht und einfach keine wirklichen Alternativen dazu.

Wer sich ernsthaft mit dem Thema Templatemetaprogrammierung auseinandersetzen will, dem würde ich einige Pflichtlektüren empfehlen. "Modern C++ Design" von Andrei Alexandrescu ist vielleicht ein guter, wenn auch nicht einfacher Einstieg. Dort werden so gut wie alle der grundlegenden Techniken erklärt, aber man sollte schon recht sicher im Umgang mit C++ sein. Das Buch "C++ Template Metaprogramming" von David Abrahams und Aleksey Gurtovoy erklärt die Thematik recht anschaulich am Beispiel von Boost.MPL, ohne allerdings detailliert die Implementierungen der einzelnen Algorithemen zu erläutern. Und zu guter Letzt sei noch das Buch "C++ Templates" von David Vandevoorde und Nicolai Josuttis erwähnt, welches sich ausführlichst mit Templates beschäftigt und damit die (theoretischen) Grundlagen für vielen Methoden liefert.

Der Artikel ist wie folgt gegliedert. Nach dieser Einführung werde ich im zweiten Teil zeigen, wie numerische Konstanten für Templatemetaprogrammierung definiert werden. Wir werden diese Definition später, wenn es um die Konstruktion der numerische Algorithmen geht, wieder benötigen. Danach werde ich anhand eines Beispieles kurz erläutern, welche Probleme generell beim Erstellen von Templatemetaprogrammen auftreten können, wie man am besten Fehler in diesen Programmen findet und Templatemetaprogramme debuggt. Danach kommt der Hauptteil, in welchem die klassischen Runge-Kutta-Verfahren zum Lösen von Differentialgleichungen mit Hilfe von Templatemetaprogrammen konstruiert werden. Der Artikel endet mit einer Diskussion der Ergebnisse und Methoden.

Numerische Konstanten

Numerische Konstanten können mit den eingebauten integralen Typen der MPL dargestellt werden:

typedef mpl::int_< 0 > zero;
typedef mpl::char_< 1 > one;
typedef mpl::long_< 2 > two;
typedef mpl::long_< -2 > minus_two;

Der numerische Wert wird gewöhnlich in einer statischen Membervariable ::value abgelegt und kann wie folgt ausgegeben werden.

cout << zero::value << endl;
cout << one::value << endl;
cout << two::value << endl;
cout << minus_two::value << endl;

Ausserdem sind diese Typen nullary Metafunktionen, welche ihren Typ zurückgeben. Daher sind die beiden Typen mpl::int_< 1 > und mpl::int_< 1 >::type gleichwertig.

Man beachte: Da Zahlen als Typen definiert sind, steht ihr Wert während des Kompiliervorgangs fest. Ausdrücke, welche diese Werte benutzen, können dann vom Kompiler optimiert werden. Allerdings können nur die ganzzahligen Typen (char, short, int, long, ...) mit Hilfe der numerischen Konstanten der MPL ausgedrückt werden. Benötigen wir Fließkommazahlen zur Compilezeit, so können wir diese durch rationale Zahlen darstellen:

template< class T , T Numerator , T Denominator >
struct rational_number
{
    typedef T value_type;
    typedef rational_number type;
    const static value_type numerator = Numerator;
    const static value_type denominator = Denominator;
};

template< class Num >
struct conv
{
    static double get_value( void )
    {
        return static_cast< double >( Num::value );
    }
};

template< class T , T Numerator , T Denominator >
struct conv< rational_number< T , Numerator , Denominator > >
{
    typedef rational_number< T , Numerator , Denominator > num_type;
    static double get_value( void )
    {
        return static_cast< double >( num_type::numerator ) / static_cast< double >( num_type::denominator );
    }
};

Die Klasse conv übernimmt hier nur die Konvertierung zu double s. Man sieht aber hier auch einen bedeutenden Unterschied zu integralen Typen. Der Wert einer rationalen Zahl kann nicht so ohne Weiteres in einer statischen Member-Variable abgelegt werden, da dies nur für integrale Typen erlaubt ist, es sein denn, man definiert diese Variable extern, was aber nicht besonders ansprechend ist. Außerdem ist die Variable dann auch wirklich als Symbol vorhanden. Daher muss hier der Umweg über eine statische Klasse gegangen werden. Diese Klasse ist spezialisiert für rational_number und man kann mit ihr sowohl integrale als auch Fließkommazahlen einheitlich benutzen, beispielsweise so:

typedef rational_number< long , 10 , 15 > two_third;
cout << conv< two_third >::get_value() << endl;
cout << conv< two >::get_value() << endl;

Für unsere Anwendung ist die Definition von numerischen Konstanten völlig ausreichend. Zusätzlich könnte man jetzt noch Operationen (+, -, *, /) definieren. Es gibt sogar einige Bibliotheken, welche die Benutzung von numerischen Konstanten soweit auf die Spitze treiben, dass beispielsweise trigonometrische Funktionen wie der Sinus oder der Cosinus vom Compiler berechnet werden können, siehe hier.

Beispielalgorithmen und Debugging

Debugging von Templatemetaalgorithmen ist meistens problematisch. Eine richtige Lösung hab ich nicht, allerdings kann man mit Hilfe von typeid einige Informationen über mögliche Probleme erhalten. Nehmen wir ein Beispiel: Wir wollen einen beliebig langen fusion -Vektor konstruieren, dessen Elemente tr1::array s sind. Die Länge der Arrays soll dabei ihrem Index in dem Vektor entsprechen. Beispielsweise soll ein Vektor der Länge 3 so aussehen

typedef fusion::vector
<
    tr1::array< double , 1 > ,
    tr1::array< double , 2 > ,
    tr1::array< double , 3 >
> coef_type;

Ein erster Versuch solch einen Vektor zu konstruieren wäre

typedef mpl::copy
<
    mpl::range_c< size_t , 1 , N + 1 > ,
    mpl::inserter
    <
        mpl::vector<> , 
        mpl::push_back
        <
            mpl::_1 ,
            tr1::array< double , mpl::_2::value >
        >
    >
>::type mpl_coef_type;
typedef fusion::result_of::as_vector< mpl_coef_type >::type coef_type;

Dies würde einem klassischen C++ Kopieralgorithmus mit einem speziellen back_insert_iterator entsprechen, ganz analog zu

// Achtung, das Beispiel wird nicht kompilieren
vector< size_t > indices;
for( size_t i=1 ; i<=N ; ++i ) indices.push_back( i );
vector< arrays > vec;
copy( indices.begin() , indices.begin() , back_insert_iterator< vector< arrays > >( vec ) );

Aber gehen wir Schritt für Schritt durch den Meta-Algorithmus:

1. mpl::copy< Seq , Inserter > ist der zentrale Kopiervorgang. Er nimmt eine MPL-Sequenz Seq entgegen und gibt jedes Element dieser Sequenz an den Inserter weiter.
2. mpl::range_c< size_t , 1 , N + 1 > konstruiert einen mpl::vector dessen Elemente size_t s von 1 bis N sind, also wie mpl::vector< 1 , 2 , 3 , ... , N-1 , N > . Dieser Vektor wird benötigt, da er erstens die richtige Anzahl von Elementen für den finalen Vektor besitzt und zweitens jedes Element die Länge des aktuellen Arrays ist.
3. mpl::inserter< Seq, Op > ist ein allgemeiner Inserter, ähnlich zu den Insert-Iteratoren der STL. Als erstes Template-Argument wird hier eine initiale Sequenz erwartet, in unserem Falle ein leerer Vektor mpl::vector<> . Das zweite Argument definiert die Einfüge-Operation.
4. mpl::push_back< Seq , El > macht nichts anderes als El an die Sequenz Seq anzuhängen. In unserem Fall wird also an den aktuellen Vektor ein neues Element vom Typ std::tr1::array angehängt. Hier sieht man auch sehr schön die Verwendung von Platzhaltern. mpl::_1 und mpl::_2 bezeichnen das erste bzw. zweite Argument, welches an push_back übergeben wird.
5. fusion::result_of::as_vector konvertiert den Mpl-Vektor in einen Fusion-Vektor.

Ausgeschrieben soll wird der gesamte Algorithmus in etwa so aussehen:

mpl::push_back
<
    mpl::push_back
    < 
        mpl::push_back
        <
            mpl::vector<> ,
            tr1::array< double , 1 >
        >::type ,
        tr1::array< double , 2 >
    >::type ,
    tr1::array< double , 3 >
>::type

Der Code kompiliert und man könnte meinen, dass wir die Aufgabe geschafft haben. Aber dem ist nicht so. Wollen wir diesen jetzt verwenden, wird der Compiler eine schier unglaubliche Anzahl von Fehlern ausspucken. Diese nach ihrer Essenz zu durchforsten, ist mühsam und selbst wenn man die entscheidenden Stellen gefunden hat, kann es sein, dass man mit der Fehlermeldung nicht viel anfangen kann (mir geht das zumindest so). Um diese Fehlersuche abzukürzen, gibt es mehrere Wege:

1. Benutze STLFilt! Es wird die Anzahl der Fehlermeldungen drastisch reduzieren.
2. Benutze statische Assertions! Damit können bestimmte Fehler viel früher gefunden werden. Statische Assertions sind insbesondere für Library-Entwickler wichtig, weil damit Bedingungen an die Typen geprüft werden.
3. Benutze einen Typdebugger! Er wird Information über den Typ liefern, welche eventuell hilfreich sind.

template< class T >
struct wrap {};

struct debug_type
{
    size_t m_indent;

    debug_type( size_t indent = 0 )
    : m_indent( indent ) { }

    template< class T >
    void operator()( wrap<T> ) const
    {
        typedef typename boost::mpl::is_sequence< T >::type is_seq;
        debug_impl( wrap< T >() , is_seq() );
    }

    template< class Seq >
    void debug_impl( wrap< Seq > , boost::mpl::true_ ) const
    {
        std::cout << indent() << "Sequence : " << typeid(Seq).name() << std::endl;
        boost::mpl::for_each<
        typename boost::mpl::transform< Seq , wrap< boost::mpl::_1 > >::type >( debug_type( m_indent + 1 ) );
    }

    template< class T >
    void debug_impl( wrap< T > , boost::mpl::false_ ) const
    {
        std::cout << indent() << typeid(T).name() << std::endl;
    }

    std::string indent( void ) const
    {
        std::string in;
        for( size_t i=0 ; i<m_indent ; ++i ) in += "  ";
        return in;
    }
};

template< class T >
void debug( size_t indent = 0 )
{
    debug_type d( indent );
    wrap< T > t;
    d( t );
}

Der wrap -Trick kommt aus dem Buch von David Abrahams und sorgt dafür, dass der Typ-Printer auch mit Referenzen funktioniert (Referenzen sind nicht default-konstruierbar und können nicht einfach von mpl::for_each erzeugt werden).

In unserem Beispiel des Vektors von Arrays führt nun die Benutzung von

debug< mpl_coef_type >();
debug< coef_type >();
dout << endl;

zu folgender Ausgabe:

Sequence : N5boost3mpl6v_itemINSt3tr15arrayIdLm2EEENS1_IS4_NS1_IS4_NS0_6vectorIN4mpl_2naES7_S7_S7_S7_S7_S7_S7_S7_S7_S7_S7_S7_S7_S7_S7_S7_S7_S7_S7_EELi0EEELi0EEELi0EEE
  NSt3tr15arrayIdLm2EEE
  NSt3tr15arrayIdLm2EEE
  NSt3tr15arrayIdLm2EEE
Sequence : N5boost6fusion7vector3INSt3tr15arrayIdLm2EEES4_S4_EE
  NSt3tr15arrayIdLm2EEE
  NSt3tr15arrayIdLm2EEE
  NSt3tr15arrayIdLm2EEE

Und voilà, man sieht, dass mpl::_2::value nicht richtig aufgelöst/ersetzt wurde, denn die Arrays haben die Länge 2. Eine Lösung ist schnell gefunden, wir packen die Definition des Arrays einfach in eine Meta-Funktion:

template< class T , class Constant >
struct array_wrapper
{
    typedef typename tr1::array< T , Constant::value > type;
};

typedef mpl::copy
<
    mpl::range_c< size_t , 1 , N + 1 > ,
    mpl::inserter
    <
        mpl::vector<> ,
        mpl::push_back
        <
            mpl::_1 ,
            array_wrapper< double , mpl::_2 >
        >
    >
>::type mpl_coef_type;
typedef fusion::result_of::as_vector< mpl_coef_type >::type coef_type;

und unser kleiner, aber feiner Algorithmus liefert das gewünschte Ergebnis.

Konstruktion von numerischen Algorithmen

Kommen wir nun zum eigentlich Thema dieses Artikels, der Konstruktion von numerischen Algorithmen mit Hilfe von Template-Metaprogramming. Wir verdeutlichen unseren Ansatz anhand der klassischen expliziten Runge-Kutta-Solvern zum Lösen von gewöhnlichen Differentialgleichungen. Ähnliche Ansätze können auch für verwandte numerische Probleme wie

• Numerische Integration / Lagrange-Integration,
• Konstruktion von Spline-Basis-Funktionen oder auch
• Finite Differenzen zum Lösen von partiellen Differentialgleichungen

genutzt werden.

Ein gewöhnliche Differentialgleichung beschreibt die (zeitliche) Entwicklung einer Größe x anhand ihrer Ableitung(en):

dx / dt = f(x,t)

x kann dabei vektoriell sein, daher aus mehreren Werten bestehen. Existiert zu einem Zeitpunkt t=t0 eine Anfangsbedingung x(t0)=x0 und ist die zeitliche Entwicklung von x gesucht, so spricht man von einem Anfangswertproblem. Numerisch werden Anfangswertprobleme normalerweise iterativ gelöst. Man bestimmt also eine aufeinanderfolgende Sequenz x(t0),x(t0+dt),x(t0+2dt),... , welche die Lösung numerisch repäsentiert. dt ist dabei die Schrittweite.

Eine einfache Klasse von solchen Algorithmen sind sogenannte explizite Runge-Kutta-Verfahren. Sie haben den großen Vorteil, dass keine Gleichungssysteme gelöst werden müssen. Schematisch werden Runge-Kutta-Solver durch Butcher-Tableaus dargestellt:

c[t]1[/t] |
c[t]2[/t] | a[t]2,1[/t]
c[t]3[/t] | a[t]3,1[/t] a[t]3,2[/t]
...
c[t]s[/t] | a[t]s,1[/t] a[t]s,2[/t] ... a[t]s,s-1[/t]
-----------------------------------
   | b[t]1[/t] b[t]2[/t] ... b[t]s-1[/t] b[t]s[/t]

s bezeichnet dabei die Ordnung des Verfahrens, was ein wichtiger Parameter ist, da er die Genauigkeit der Lösung charakterisiert.

Aus solch einer Tabelle kann man nun das Verfahren zum Lösen der DGL wie folgt konstruieren:

k[t]1[/t] = f( x(t) , t + c[t]1[/t] * dt )
k[t]2[/t] = f( x(t) + dt * a[t]2,1[/t] * k[t]1[/t] , t + c[t]2[/t] * dt )
k[t]3[/t] = f( x(t) + dt * ( a[t]3,1[/t] * k[t]1[/t] + a[t]3,2[/t] * k[t]2[/t] ), t + c[t]3[/t] * dt )
...
k[t]s[/t] = f( x(t) + dt * ( a[t]s,1[/t] * k[t]1[/t] + ... + a[t]s,s-1[/t] * k[t]s-1[/t] ) , t + c[t]s[/t] * dt )
x(t+dt) = x(t) + dt * ( b[t]1[/t] k[t]1[/t] + ... + b[t]s[/t] k[t]s[/t] )

Ein einfaches Beispiel ist das Euler-Verfahren:

0 |
--|---
  | 1

Oder ausgeschrieben:

k[t]1[/t] = f( x(t) , t )
x(t+dt) = x(t) + dt * k[t]1[/t]

Ein etwas komplizierteres Beispiel ist das klassische Runge-Kutta-Verfahren vierter Ordnung (im Folgenden als RK4 bezeichnet):

0   |
1/2 | 1/2
1/2 | 0   1/2
1   | 0   0   1
----|----------------
    | 1/3 1/6 1/6 1/3

In Pseudocode sieht diese Methode so aus:

k[t]1[/t] = f( x(t) , t )
k[t]2[/t] = f( x(t) + dt * 1/2 * k[t]1[/t] , t + 1/2 * dt )
k[t]3[/t] = f( x(t) + dt * 1/2 * k[t]2[/t] , t + 1/2 * dt )
k[t]4[/t] = f( x(t) + dt * 1 * k[t]3[/t] , t + 1 * dt )
x(t+dt) = x(t) + dt * ( 1/6 k[t]1[/t] + 1/3 k[t]2[/t] + 1/3 k[t]3[/t] + 1/6 k[t]4[/t] )

Wir werden nun Schritt für Schritt einen Meta-Algorithmus konstruieren, welcher solche Algorithmen während des Kompiliervorganges erzeugt. Als erstes werden die Koeffizienten a[t]i,j[/t] , b[t]i[/t] und c[t]i[/t] definiert. Danach wird das Hauptgerüst des Algorithmus entworfen und im letzten Teil die einzelnen Schritte ausgeführt. Zum Schluss wird eine alternative Implementation vorgestellt, welche sich durch eine einfache Konstrukion der Koeffizienten auszeichnet.

Die Koeffizienten eines Verfahrens werden in MPL-Vektoren abgelegt. Für RK4 sieht die Definition wie folgt aus:

// allgemeine Konstanten
typedef mpl::int_< 0 > null;
typedef mpl::int_< 1 > one;
typedef rational_number< long , 1 , 2 > one_half;
typedef rational_number< long , 1 , 3 > one_third;
typedef rational_number< long , 1 , 6 > one_sixth;

typedef mpl::vector
<
    mpl::vector< one_half > ,
    mpl::vector< null , one_half > ,
    mpl::vector< null , null , one >
> rk4_a;
typedef mpl::vector< null , one_half , one_half , one > rk4_c;
typedef mpl::vector< one_sixth , one_third , one_third , one_sixth > rk4_b;

Der Solver selbst ist eine templatisierte Klasse, welche durch die Koeffizienten parametrisiert ist:

template< class State , class ButcherA , class ButcherB , class ButcherC >
class explicit_runge_kutta
{
public:
    typedef ButcherA butcher_a;
    typedef ButcherB butcher_b;
    typedef ButcherC butcher_c;

    static const size_t num_of_sages = mpl::size< ButcherC >::value;

    typedef double time_type;
    typedef State state_type;
...
}

time_type und state_type geben die Typen von t bzw. x an und für state_type wird ein std::tr1::array< double , N > erwartet.

Als Nächstes überprüfen wir, ob butcher_a , butcher_b und butcher_c den Anforderungen an die Koeffizienten genügen. Insbesondere prüfen wir, ob diese Typlisten die entsprechenden Anzahlen von Elementen enthalten.

template< class State , class ButcherA , class ButcherB , class ButcherC >
class explicit_runge_kutta
{
    ...
    BOOST_STATIC_ASSERT( ( num_of_stages > 0 ) );
    BOOST_STATIC_ASSERT( ( num_of_stages == size_t( mpl::size< butcher_b >::value ) ) );
    BOOST_STATIC_ASSERT( ( ( num_of_stages - 1 ) == size_t( mpl::size< butcher_a >::value ) ) ) ;

    typedef typename mpl::equal
    <
        mpl::range_c< int , 1 , num_of_stages > ,
        typename mpl::transform< butcher_a , mpl::size< mpl::_1 > >::type ,
        mpl::equal_to< mpl::_1 , mpl::_2 >
    >::type equal_type;
    BOOST_STATIC_ASSERT(( equal_type::value == true ));
};

Jetzt kommen wir zum eigentlichen Kernstück unseres Metaalgorithmus. Die expliziten Runge-Kutta-Verfahren verfügen über s Stufen, wobei s hier wieder die Ordnung des Verfahrens ist. Die erste Stufe sieht wie folgt aus:

k[t]1[/t] = f( x(t) , t + c[t]1[/t] )
xtmp = x(t) + dt * a[t]2,1[/t] * k[t]1[/t]

In den nachfolgenden Stufen (außer der letzten) wird immer

k[t]i[/t] = f( xtmp , t + c[t]i[/t] )
xtmp = x(t) + dt * a[t]i+1,1[/t] * k[t]1[/t] + ... + a[t]i+1,i[/t] * k[t]i[/t]

berechnet. Die letzte Stufe unterscheidet sich von den anderen, hier wird das Ergebnis geschrieben:

k[t]s[/t] = f( xtmp , t + c[t]s[/t] )
x(t+dt) = x(t) + dt * b[t]1[/t] * k[t]1[/t] + ... + b[t]s[/t] * k[t]s[/t]

Diesen Algorithmus setzten wir nun in einen Meta-Algorithmus um:

template< class State , class ButcherA , class ButcherB , class ButcherC >
class explicit_runge_kutta_mpl
{
public:
    ...
    typedef mpl::range_c< size_t , 0 , num_of_stages > indices;
    typedef typename mpl::push_back< butcher_a , butcher_b >::type butcher_right;
    typedef mpl::zip_view< mpl::vector< indices , butcher_c , butcher_right > > butcher_tableau;

    template< class System >
    void do_step( System &system , state_type &x , double t , const double dt )
    {
        mpl::for_each< butcher_tableau >( calculate_stage< System >( system , x , m_x_tmp , m_k_vector , t , dt ) );
    }
private:

    state_type m_k_vector[num_of_stages];
    state_type m_x_tmp;
    ...
};

Schauen wir uns an, was hier passiert. butcher_tableau ist ein Mpl-Vektor mit folgenden Einträgen:

0 , c[t]1[/t] , a[t]2,1[/t] 
1 , c[t]2[/t] , a[t]3,1[/t] , a[t]3,2[/t]
2 , c[t]3[/t] , a[t]4,1[/t] , a[t]4,2[/t] , a[t]4,3[/t] 
...
s , c[t]s[/t] , b[t]1[/t]   , b[t]2[/t]   , b[t]3[/t]   , ... , b[t]s[/t]

Die i-te Zeile enthält also alle Informationen, welche wir zur Berechnung der i-ten Stufe benötigen. Über diesen Vektor wird nun in der Zeile

mpl::for_each< butcher_tableau >( calculate_stage< System >( system , x , m_x_tmp , m_k_vector , t , dt ) );

iteriert und in jeder Iteration calculate_state aufgerufen. Die Methode do_step geht nun einen Zeitschritt in der Differentialgleichung nach vorn, do_step berechnet x(t+dt) anhand von x(t) . Die Differentialgleichung selbst wird durch system angegeben. system kann hier ein Funktor, aber auch ein Funktionspointer sein. Es wird dabei immer erwartet, dass wir

state_type dxdt;
system( x , dxdt , t );

ausführen können.

Nun wenden wir uns der der eigentlichen Berechnung jeder einzelnen Stufe zu. Das passiert in calculate_stage , welche eine Member-Struktur von explicit_runge_kutta_mpl ist. Diese Struktur ist gleichzeitig auch ein Funktor, daher existiert der Klammer-Operator.

template< class State , class ButcherA , class ButcherB , class ButcherC >
class explicit_runge_kutta_mpl
{
public:
    ...
    template< class System >
    struct calculate_stage
    {
        System &system;
        state_type &x , &x_tmp;
        state_type *k_vector;
        const double t;
        const double dt;

        calculate_stage( System &_system , state_type &_x , state_type &_x_tmp , state_type *_k_vector , const double _t , const double _dt )
        : system( _system ) , x( _x ) , x_tmp( _x_tmp ) , k_vector( _k_vector ) , t( _t ) , dt( _dt )
        {}

        template< class Stage >
        void operator()( Stage v )
        {
            typedef typename mpl::at< Stage , mpl::int_< 0 > >::type index_type;
            typedef typename mpl::at< Stage , mpl::int_< 1 > >::type time_value_type;
            typedef typename mpl::at< Stage , mpl::int_< 2 > >::type coef_type;

            const static size_t index = index_type::value;
            const static double time_value = conv< time_value_type >::get_value();

            BOOST_STATIC_ASSERT(( ( mpl::size< coef_type >::value - 1 ) == index ));

            if( index == 0 ) system( x , k_vector[index] , t + time_value * dt );
            else system( x_tmp , k_vector[index] , t + time_value * dt );

            if( index == ( num_of_stages - 1 ) )
                algebra< state_type , coef_type , mpl::int_< index > >::do_step( x , x , k_vector , dt );
            else
                algebra< state_type , coef_type , mpl::int_< index > >::do_step( x_tmp , x , k_vector , dt );
        }
    };
    ...
};

Das Klassen-Template algebra enthält nur eine einzelne statische Methode do_step und ist durch den Zustandstyp state_type , die Koeffizienten des Verfahrens coef_type und die Ordnung des Verfahrens order parametrisiert.

template< class state_type , class coef_tye , class order > struct algebra_mpl;

algebra muss für alle Ordnungen spezialisiert werden. Hier könnte bestimmt der Präprozessor helfen, aber das ist eine andere Geschichte. Deswegen schreiben wir alle Ordnungen aus. Zum Beispiel sehen die ersten beiden Ordnungen so aus:

template< class state_type , class coef_type >
struct algebra_mpl< state_type , coef_type , mpl::int_< 0 > >
{
    static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , const double dt )
    {
        const static double a1 = dt * conv< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value();

        for( size_t i=0 ; i<x.size() ; ++i )
        {
            x_tmp[i] = x[i] + a1 * k_vector[0][i];
        }
    }
};

template< class state_type , class coef_type >
struct algebra_mpl< state_type , coef_type , mpl::int_< 1 > >
{
    static void do_step( state_type &x_tmp , const state_type &x , const state_type *k_vector , const double dt )
    {
        const static double a1 = dt * conv< typename mpl::at< coef_type , mpl::int_< 0 > >::type >::get_value();
        const static double a2 = dt * conv< typename mpl::at< coef_type , mpl::int_< 1 > >::type >::get_value();

        for( size_t i=0 ; i<x.size() ; ++i )
        {
            x_tmp[i] = x[i] + a1 * k_vector[0][i] + a2 * k_vector[1][i];
        }
    }
};

Damit sind wir fertig. Das Letzte, was uns jetzt noch übrig bleibt, ist einen speziellen Solver zu initialisieren sowie zu instanzieren und dann können wir mit dem Lösen von Differentialgleichungen beginnen. Für das vorangegangene Euler-Beispiel sieht die Definition des Solvers wie folgt aus:

typedef mpl::int_< 0 > null;
typedef mpl::int_< 1 > one;

typedef mpl::vector< > euler_a;
typedef mpl::vector< one > euler_b;
typedef mpl::vector< null > euler_c;

template< class State > class stepper_euler
: public explicit_runge_kutta_mpl< State , euler_a , euler_b , euler_c > {};

Als erstes werden die Konstanten für das Verfahren definiert, in diesem Fall 0 und 1. Danach werden die Koeffizienten a[t]i,j[/t] , b[t]i[/t] und c[t]i[/t] definiert. Die Koeffizienten b und c enthalten jeweils nur einen Eintrag, die Koeffizienten a sind sogar leer. Danach wird der Euler-Solver per Vererbung definiert. Für den klassischen Runge-Kutta-Solver sieht alles sehr ähnlich aus, nur dass andere Koeffizienten vorkommen:

typedef rational_number< int , 1 , 2 > one_half;
typedef rational_number< int , 1 , 3 > one_third;
typedef rational_number< int , 1 , 6 > one_sixth;

typedef mpl::vector
<
    mpl::vector< one_half > ,
    mpl::vector< null , one_half > ,
    mpl::vector< null , null , one >
> rk4_a;
typedef mpl::vector< null , one_half , one_half , one > rk4_c;
typedef mpl::vector< one_sixth , one_third , one_third , one_sixth > rk4_b;

template< class State > class stepper_rk4
: public explicit_runge_kutta_mpl< State , rk4_a , rk4_b , rk4_c > { };

Alternative Implementierung

Ein Nachteil der obigen Implementierung des Compilezeit-Runge-Kutta Algorithmus ist sicherlich die Art und Weise, wie die Koeffizienten angegeben werden müssen. Für einen Koeffizienten wie 1/6 müssen wir beispielsweise rational_number< int , 1 , 6 > schreiben. Eine einfachere Definition könnte beispielsweise konstante Arrays benutzen, welche im Konstruktor des Runge-Kutta-Verfahrens übergeben werden. Wir werden diesen Ansatz jetzt entwickeln, er erfordert allerdings ein kleines bisschen mehr Metaprogrammierungszauberei.

In der alternativen Implementierung werden die Typen der Koeffizienten folgendermaßen definiert:

template< class T , class Constant >
struct array_wrapper
{
    typedef typename std::tr1::array< T , Constant::value > type;
};

template< class StateType , size_t NumOfStages >
class explicit_runge_kutta_fusion
{
public:
    typedef StateType state_type;
    const static size_t num_of_stages = NumOfStages;

    typedef typename fusion::result_of::as_vector
    <
        typename mpl::copy
        <
            mpl::range_c< size_t , 1 , num_of_stages > ,
            mpl::inserter
            <
                mpl::vector0< > ,
                mpl::push_back< mpl::_1 , array_wrapper< double , mpl::_2 > >
            >
        >::type
    >::type coef_a_type;

    typedef std::tr1::array< double , num_of_stages > coef_b_type;
    typedef std::tr1::array< double , num_of_stages > coef_c_type;
    ...
}

Die Koeffizienten b und c werden in einfachen Arrays gespeichert. Für die Koeffizienten a benötigen wir allerdings etwas Komplizierteres: einen fusion::vector . Solche Vektoren sind ähnlich zu std::tr1::tuple , aber mit dem Unterschied, dass dazu eine riesige Menge an Meta-Algorithmen und Funktionen zur Verfügung stehen. Der Meta-Algorithmus hier konstruiert einen solchen Vektor (siehe auch das obige Beispiel). result_of::as_vector konvertiert einen MPL-Vektor in einen Fusion-Vektor. Die Einträge in diesem Vektor sind auch hier std::tr1::array s, wobei die Länge der Arrays mit dem Index im Vektor zunimmt, daher sieht der Vektor wie folgt aus:

fusion::vector< std::tr1::array< double , 1 > , std::tr1::array< double , 2 > , std::tr1::array< double , 3 > , ... >

So vereinfacht sich dann beispielsweise die Koeffizienten-Definition des klassischen Runge-Kutta-Verfahrens zu

using namespace std::tr1;
const array< double , 1 > coef_a1_rk4 = {{ 0.5 }};
const array< double , 2 > coef_a2_rk4 = {{ 0.0 , 0.5 }};
const array< double , 3 > coef_a3_rk4 = {{ 0.0 , 0.0 , 1.0 }};
const fusion::vector< array< double , 1 > , array< double , 2 > , array< double , 3 > > coef_a_rk4( coef_a1_rk4 , coef_a2_rk4 , coef_a3_rk4 );
const array< double , 4 > coef_b_rk4 = {{ 1.0 / 6.0 , 1.0 / 3.0 , 1.0 / 3.0 , 1.0 / 6.0 }};
const array< double , 4 > coef_c_rk4 = {{ 0.0 , 0.5 , 0.5 , 1.0 }};

Auch hier wird der eigentliche Algorithmus in mehrere Stages unterteilt. Diese sind exakt äquivalent zur MPL-Implementierung, außer dass die Koeffizienten nicht als Typen gespeichert werden, sondern als Array bzw. in Fusion-Vektoren. Der entscheidende Punkt ist dabei, dass diese Arrays/Vektoren als const deklariert werden und die Längen fix sind. Dadurch kann der Compiler (hoffentlich) Optimierungen anwenden, so als ob der Algorithmus direkt ausgeschrieben würde.

Der Rest der Implementierung ist dann:

template< class T , class Constant >
struct stage_fusion_wrapper
{
    typedef typename fusion::vector< size_t , T , std::tr1::array< T , Constant::value > > type;
};

template< class StateType , size_t NumOfStages >
class explicit_runge_kutta_fusion
{
    ...

    typedef typename fusion::result_of::as_vector
    <
        typename mpl::copy
        <
            mpl::range_c< int , 1 , num_of_stages + 1 > ,
            mpl::inserter
            <
                mpl::vector0<> ,
                mpl::push_back< mpl::_1 , stage_fusion_wrapper< double , mpl::_2 > >
            >
        >::type
    >::type stage_vector_base;

    struct stage_vector : public stage_vector_base
    {
        struct do_insertion
        {
            stage_vector_base &m_base;
            const coef_a_type &m_a;
            const coef_c_type &m_c;

            do_insertion( stage_vector_base &base , const coef_a_type &a , const coef_c_type &c )
            : m_base( base ) , m_a( a ) , m_c( c ) { }

            template< class Index >
            void operator()( Index ) const
            {
                fusion::at_c< 0 >( fusion::at< Index >( m_base ) ) = Index::value;
                fusion::at_c< 1 >( fusion::at< Index >( m_base ) ) = m_c[ Index::value ];
                fusion::at_c< 2 >( fusion::at< Index >( m_base ) ) = fusion::at< Index >( m_a );
            }
        };

        stage_vector( const coef_a_type &a , const coef_b_type &b , const coef_c_type &c )
        {
            typedef mpl::range_c< size_t , 0 , num_of_stages - 1 > indices;
            mpl::for_each< indices >( do_insertion( *this , a , c ) );
            fusion::at_c< 0 >( fusion::at_c< num_of_stages - 1 >( *this ) ) = num_of_stages - 1 ;
            fusion::at_c< 1 >( fusion::at_c< num_of_stages - 1 >( *this ) ) = c[ num_of_stages - 1 ];
            fusion::at_c< 2 >( fusion::at_c< num_of_stages - 1 >( *this ) ) = b;
        }
    };

    template< class System >
    struct calculate_stage
    {
        System &system;
        state_type &x , &x_tmp;
        state_type *k_vector;
        const double t;
        const double dt;

        calculate_stage( System &_system , state_type &_x , state_type &_x_tmp , state_type *_k_vector , const double _t , const double _dt )
        : system( _system ) , x( _x ) , x_tmp( _x_tmp ) , k_vector( _k_vector ) , t( _t ) , dt( _dt ) { }

        template< typename T , size_t stage_number >
        inline void operator()( fusion::vector< size_t , T , std::tr1::array< T , stage_number > > const &stage ) const
        {
            double c = fusion::at_c< 1 >( stage );
            if( stage_number == 1 )
                system( x , k_vector[stage_number-1] , t + c * dt );
            else
                system( x_tmp , k_vector[stage_number-1] , t + c * dt );

            if( stage_number == num_of_stages )
                fusion_algebra<stage_number>::foreach( x , x , fusion::at_c< 2 >( stage ) , k_vector , dt);
            else
                fusion_algebra<stage_number>::foreach( x_tmp , x , fusion::at_c< 2 >( stage ) , k_vector , dt);
        }
    };

public:

    explicit_runge_kutta_fusion( const coef_a_type &a , const coef_b_type &b , const coef_c_type &c )
    : m_a( a ) , m_b( b ) , m_c( c ) , m_stages( a , b , c ) , m_x_tmp() , m_k_vector() { }

    template< class System >
    inline void do_step( System &system , state_type &x , double t , const double dt )
    {
        fusion::for_each( m_stages , calculate_stage< System >( system , x , m_x_tmp , m_k_vector , t , dt ) );
    }

private:

    const coef_a_type m_a;
    const coef_b_type m_b;
    const coef_c_type m_c;
    const stage_vector m_stages;
    state_type m_x_tmp;
    state_type m_k_vector[num_of_stages];
};

Die einzelnen Stages werden hier auch wieder in einer eigenen Struktur gespeichert: stage_vector . Der Funktor calculate_stage übernimmt die Berechnung einer einzelnen Stage.

Performance

Die Performance der beiden Algorithmen ist etwa äquivalent. Außerdem stehen beide Algorithmen den ausgeschriebenen Varianten in nichts hinterher, in einigen Fällen war der Fusion-Algorithmus sogar schneller als der ausgeschriebene. Genaue Werte sind für verschiedene Compiler unterschiedlich, es gibt sogar Unterschiede zwischen gcc 4.4 und gcc 4.5.

Diskussion

Die Vorteile der beiden Ansätze liegen sofort auf der Hand. Man hat nur genau einen Algorithmus zu warten. Damit können praktisch alle klassischen expliziten Runge-Kutta-Solver (wie Euler, RK4, Midpoint, Dormand-Prince 5-4, Dormand Prince 8-5-3, Fehlberg, ...) geschrieben werden. Außerdem stehen diese Algorithmen den ausgeschriebenen Algorithmen in Performance in nichts nach. Man kann dieses Konzept natürlich auch einfach auf die Stepper mit Fehlerberechnung oder auf Stepper mit kontinuierlichem Output erweitern.

Die Nachteile bei dieser Art von Programmierung sind gering. Der Code ist schwerer zu lesen als klassischer C/C++-Code und es wird auch relativ viel Boilerplate erzeugt. Allerdings wird die dadurch längere Entwicklungszeit durch die Allgemeinheit der Verfahren mehr als wettgemacht.

Selbstverständlich erheben die beiden Versionen des Templatemetaprogrammierungs-Runge-Kutta-Verfahrens keinen Anspruch auf Allgemeingültigkeit. Sie wurden im Rahmen von Boost.Odeint entwickelt und werden dort zum Einsatz kommen. Die Erweiterungsmöglichkeiten dazu sind vielzählig:

• Abstrahierung des Zustandstyps. In den obigen beiden Beispielen wird der Zustand der Differentialgleichung durch std::tr1::array< double , N > beschrieben. Das kann selbstverständlich verallgemeinert werden.
• Benutzungen von beliebigen Floating-Point-Typen wie float , std::complex< double > , GMP, ...
• Erweiterung auf Solver mit Fehlerabschätzung und Dense-Output für Schrittweitensteuerung.
• Verallgemeinerung der Algebra.
• Höhere Ordnungen der Algebra. Im Moment werden nur Stepper bis Ordnung 4 unterstützt. Das ist trivial erweiterbar auf höhere Ordnungen. Auch können diese Erweiterungen durch Präprozessormakros automatisch erzeugt werden.

Über Kritik, Anregungen, Hinweise sowie Verbesserungsvorschläge würde ich mich sehr freuen.

Source code

Der Source code zu diesem Artikel kann hier gefunden werden.

Literatur

David Abrahams, Aleksey Gurtovoy, C++ Template Metaprogramming: Concepts, Tools, and Techniques from Boost and Beyond (Addison-Wesley, 2004).

Andrei Alexandrescu, Modern C++ Design, Generic Programming and Design Patterns Applied (Addison-Wesley, 2001)

Press William H et al., Numerical Recipes 3rd Edition: The Art of Scientific Computing, 3rd ed. (Cambridge University Press, 2007).

Ernst Hairer, Syvert P. Nørsett, and Gerhard Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems, 2nd ed. (Springer, Berlin, 2009).

Milton Abramowitz, Irene A. Stegun, Handbook of Mathematical Functions (New York: Dover Publications, 1972). (online verfügbar)

[GPC]

Ich hab's zwar schon intern in der Redaktion gesagt, aber ich sag's gern nochmal: Sehr gute Arbeit!

[Tobias Gerg]

Von mir auch ein großes Lob... Eine sehr guter Artikel...

[mariomulansky]

kleine anmerkung noch dazu:
als erstes hatten wir eine Implementation des generischen Runge-Kuttas, die komplett "dynamisch" umgesetzt war, also ohne Metaprogramming funktionierte. Da war das Butcher-Tableau ein vector< vector < double > > und es wurde in 2 for-schleifen ueber die koeffizienten iteriert.
Das war nur leider deutlich (faktor 3) langsamer als eine direkte Implementation und daher die Template-Version.
Das vielleicht noch als weitere Motivation warum hier Metaprogramming benutzt wurde.

[knivil]

Vielleicht ganz interessant und zum Thema passend: Haskell and C++ Template Metaprogramming. Zwar ist der Vortragende seeehr langsam und kritisiert C++ Template Metaprogramming oft genug, aber nicht umsonst sind andere wissenschaftliche Bibliotheken nicht in C++ geschrieben. Performance ist mittlerweile zweitrangig, Parallelisierbarkeit steht an erster Stelle. Leider ist das bei Differentialgleichungen nicht immer gegeben. Nichtsdestotrotz ist es ein grossartiger Artikel.

Zu C++ Concepts und Haskell gibt es auch einige Beitraege auf der verlinkten Seite. Wer mag kann auch gern die type class SearchProblem aus Escape from Zurg als Template zur Uebung formulieren.

[headmyshoulder]

Ich hab mir gerade das Video angesehen, zumindest zum Teil, und Haskell sieht sehr viel einfacher aus. Die Art und Weise wie Haskel smit Spezialisierungen umgeht, erscheint mit sehr elegant aus und macht den Code sehr schlank. Obwohl er auch schon einige Negativbeispiel der C++ Welt gezeigt hat. Der Anwendungsbereich von TMP ist halt auch nicht so riesig, und meist geht es dabei um Performance oder Typprüfungen.

Das Beispiel der Fakultät kann wahrscheinlich jeder gute Compiler mit der Standardimplementierung

int fac( int N ) { return ( N == 0 ) ? 1 : fac( N - 1 );

lösen. Falls N als Konstante übergeben wird, ersetzt der Compiler das mit der Lösung.

[314159265358979]

Ich finde den Artikel sehr gut, aber wenn man von den Themen noch nie vorher gehört hat wie ich, ist das Thema zu komplex, da man sich ziemlich lange einlesen muss.

[abc.w]

Hallo, ich habe eine Frage. Sehe ich das richtig, dass die beschriebene Numerik und Template-Metaprogrammierung nur für die Eingabedaten funktionieren, die dem Compiler beim Bauen bekannt sind... und was ist dann der praktische Sinn dieser Methode oder dieser Art der Programmierung Man hat diesen einen schönen Algorithmus, aber er funktioniert nur beim Bauen oder wie

[headmyshoulder]

abc.w schrieb:

Hallo, ich habe eine Frage. Sehe ich das richtig, dass die beschriebene Numerik und Template-Metaprogrammierung nur für die Eingabedaten funktionieren, die dem Compiler beim Bauen bekannt sind... und was ist dann der praktische Sinn dieser Methode oder dieser Art der Programmierung Man hat diesen einen schönen Algorithmus, aber er funktioniert nur beim Bauen oder wie

Auf Deine Frage: Ja und Nein . Der genaue Algorithmus, in diesem Falle die Koeffizienten, müssen dem Compiler während des Kompilierens bekannt sein. Trotzdem kann man diesen Algorithmen benutzen, auch mit Eingabedaten. Er wird die Differentialgleichung lösen, also genau das machen was Du ihm während des Kompilieren gesagt hast. Das funktioniert natürlich auch, wenn der Anfangszustand der DGL während der Laufzeit erst bestimmt wird.

Den Ansatz den ich hier vorgestellt habe, ist dafür gedacht bestimmte Algorithmen zu erzeugen. Davon gibt es nicht so viele. Mir fallen ca. 10 ein, die sinnvoll sind. Und der eigentliche Punkt ist: es soll schnell sein. Der Compiler soll und muss den genauen Algorithmus optimieren können. Deswegen kommen Ansätze mit virtuellen Funktionen oder mit variabler Länge (variabler Anzahl an Koeffizienten) nicht in Frage.