5. [R] Numerik und Template-Metaprogrammierung

[R] Numerik und Template-Metaprogrammierung

  • G

    headmyshoulder schrieb:

    GPC schrieb:

    Freu mich schon auf den Artikel, ist ja immer wieder ein spaßiges Thema, was man mit C++ Templates alles anstellen kann 🙂

    Hoffentlich kann ich Deine Erwartungen erfüllen. Ich werd mir größte Mühe geben.

    Mach dir keinen Stress 😉 Es soll ja auch Spaß machen 🙂

    0
  • G

    Hui, hier hat sich ja schon einiges getan 👍 Gib bescheid wenn man mal drüberschauen soll :xmas1:

    0
  • H

    Ja klar, mach ich. Danke. Wird aber noch ein kleines bisschen dauern. Einen Abschnitt wird es noch geben und Form, Grammatik und Rechtschreibung sind noch unter aller Sau.

    0
  • H

    Ok, ich glaub ich bin dann soweit durch mit dem Artikel. Es wäre sehr nett, wenn das jemand darüber lesen könnte, und mir vielleicht ein bisschen Feedback geben könnte. Ich bin mir nicht sicher, ob das nicht vielleicht zu schwer und unverständlich ist.

    Vielen Dank!

    0
  • G

    Wow, krasses Zeug! Das ist schon beeindruckend/beängstigend :p was man mit Template Metaprogramierung alles machen kann. Mir gefällt der Artikel gut, da er einen guten Einblick in die Materie gewährt 👍 Der Artikel ist durchaus anspruchsvoll, aber das liegt hauptsächlich an der Thematik. Verbesserungsvorschläge von mir wären:
    - Ein kleiner Abschnitt ganz am Anfang, der beschreibt wie der Artikel aufgebaut ist und warum es wichtig ist, dass z.B. das Kapitel "Numerik" ganz am Anfang kommt
    - Vielleicht ein paar Links zu den vorgestellten Algorithmen in die Wikipedia oder eine Seite deiner Wahl, damit man als nicht Mathe-Crack leichter Informationen nachholen kann
    - Abkürzungen, wie z.B. DGL für Differentialgleichung folgendermaßen einführen: "[...] numerische Verfahren zum Lösen von Differentialgleichungen (DGL) gehen."

    Ein paar kleine Rechtschreib/Satzbauprobleme werden dann beim Korrekturlesen ausgebessert, sind aber recht wenig und auch keine Kracher dabei 😉
    Ich werde noch versuchen, ein paar andere Leute zu organisieren, die ebenfalls drüberlesen 🙂

    0
  • H

    GPC schrieb:

    Wow, krasses Zeug! Das ist schon beeindruckend/beängstigend :p was man mit Template Metaprogramierung alles machen kann. Mir gefällt der Artikel gut, da er einen guten Einblick in die Materie gewährt 👍 Der Artikel ist durchaus anspruchsvoll, aber das liegt hauptsächlich an der Thematik. Verbesserungsvorschläge von mir wären:
    - Ein kleiner Abschnitt ganz am Anfang, der beschreibt wie der Artikel aufgebaut ist und warum es wichtig ist, dass z.B. das Kapitel "Numerik" ganz am Anfang kommt
    - Vielleicht ein paar Links zu den vorgestellten Algorithmen in die Wikipedia oder eine Seite deiner Wahl, damit man als nicht Mathe-Crack leichter Informationen nachholen kann
    - Abkürzungen, wie z.B. DGL für Differentialgleichung folgendermaßen einführen: "[...] numerische Verfahren zum Lösen von Differentialgleichungen (DGL) gehen."

    Ein paar kleine Rechtschreib/Satzbauprobleme werden dann beim Korrekturlesen ausgebessert, sind aber recht wenig und auch keine Kracher dabei 😉
    Ich werde noch versuchen, ein paar andere Leute zu organisieren, die ebenfalls drüberlesen 🙂

    Danke für das positive Feedback. Die Vorschläge werde ich umsetzen. Ich würde damit aber warten bis noch mehr Feedback kommt und dann alles in einem Rutsch machen. Wenn das ok ist, könnte auch einige Leute fragen die diesen Artikel reviewen, aber das ist dann vielleicht nicht so objektiv. Hier im Forum sind auch Einige welche da eventuell drüber lesen könnten.

    Ich habe ja ein kleines bisschen Angst dass das vielleicht zuviel für den Einstieg ist. Vielleicht sollte ich das noch explizit erwähnen...

    0
  • G

    headmyshoulder schrieb:

    Danke für das positive Feedback. Die Vorschläge werde ich umsetzen. Ich würde damit aber warten bis noch mehr Feedback kommt und dann alles in einem Rutsch machen. Wenn das ok ist, könnte auch einige Leute fragen die diesen Artikel reviewen, aber das ist dann vielleicht nicht so objektiv. Hier im Forum sind auch Einige welche da eventuell drüber lesen könnten.

    Kannst du gerne machen, je mehr Feedback desto besser!

    Ich habe ja ein kleines bisschen Angst dass das vielleicht zuviel für den Einstieg ist. Vielleicht sollte ich das noch explizit erwähnen...

    Wenn du sinnvolle Möglichkeiten zur Vereinfachung findest bzw. wie man sich an manchen Stellen leichter rantasten kann, dann kein Ding. Simplicity never hurt 😃

    Ich bin jetzt erstmal 'ne Woche weg, d.h. du kannst dir ruhig Zeit mit allem lassen 🙂

    0
  • P

    Toller Artikel. Hat mir sehr gut gefallen, besonders das man mal TMP in der Praxis sieht.

    Eigentlich tust du da nur lookuptabellen generieren, oder? 😉

    0
  • H

    phlox81 schrieb:

    Toller Artikel. Hat mir sehr gut gefallen, besonders das man mal TMP in der Praxis sieht.

    Danke für die Blumen. Diese Art von Programmierung lässt sich auf eine Menge andere numerischer Probleme anwenden. Immer dann wenn Ordnungen ins Spiel kommen, z.B. Splines, Numerische Integration, PDEs ... Oder auch Zustandsmachinen- und automaten, aber das ist ein anderes Thema :).

    phlox81 schrieb:

    Eigentlich tust du da nur lookuptabellen generieren, oder? 😉

    Jepp. Ausserdem werden noch die Methoden, welche einen Eintrag in der Tabelle auswerten, generiert.

    0
  • H

    Wie geht es denn jetzt weiter? Ich hab ein paar Bekannte nach Reviews gefragt, aber das kann eventuell noch ein bisschen dauern und die Objektivität ist vielleicht nicht gegeben.

    Wie lange dauert ein Review normalerweise? Könnte ich noch etwas tun um das zu beschleunigen?

    0
  • G

    headmyshoulder schrieb:

    Wie geht es denn jetzt weiter? Ich hab ein paar Bekannte nach Reviews gefragt, aber das kann eventuell noch ein bisschen dauern und die Objektivität ist vielleicht nicht gegeben.

    Wie lange dauert ein Review normalerweise? Könnte ich noch etwas tun um das zu beschleunigen?

    Du kannst gerne deine Bekannten fragen und deren Anmerkungen auch einbauen 😉 Normalerweise bestimmt der Autor, wann die "T"-Phase vorbei ist... in deinem Fall würde es sich dann anbieten, auf "R" zu wechseln, da der Artikel doch schon einige Male aufgerufen wurde, aber keine weiteren Kommentare folgten. Aber kein Stress, du kannst ruhig noch die Reviews deiner Bekannten abwarten 😉

    0
  • H

    Super, dann mach ich das so. Werd alle Kommentare einbauen und den Artikel dann auf R setzen. Danke

    0
  • H

    So, ich hab den Artikel jetzt auf R gesetzt. Ich kann natürlich trotzdem "fachliche" Anmerkungen einfügen. Ich wollte auch die Sourcen hochladen, aber irgednwie krieg ich keine Verbindung auf den ftp Server. Vielleicht sind Username/PW nicht mehr aktuell oder der ftp ist gerade down.

    0
  • G

    headmyshoulder schrieb:

    So, ich hab den Artikel jetzt auf R gesetzt. Ich kann natürlich trotzdem "fachliche" Anmerkungen einfügen. Ich wollte auch die Sourcen hochladen, aber irgednwie krieg ich keine Verbindung auf den ftp Server. Vielleicht sind Username/PW nicht mehr aktuell oder der ftp ist gerade down.

    Also ich habe gerade eben erfolgreich eine Verbindung auf den server hinbekommen. Vermutlich hast du dich vertippt 😉

    Und den Artikel werd ich heute oder morgen noch korrekturlesen. Denk auch dran, dass du noch ein Autorenprofil brauchst, in dem du ein wenig über dich schreibst. Muss ja nix großes sein 🙂

    0
  • H

    GPC schrieb:

    Also ich habe gerade eben erfolgreich eine Verbindung auf den server hinbekommen. Vermutlich hast du dich vertippt 😉

    Nimmst Du die Daten von hier? Und wenn ja, normales FTP ohne verschlüsselung? Und welcher Port? Vielleicht ist mein Client auch nicht so intelligent (filezilla).

    GPC schrieb:

    Und den Artikel werd ich heute oder morgen noch korrekturlesen. Denk auch dran, dass du noch ein Autorenprofil brauchst, in dem du ein wenig über dich schreibst. Muss ja nix großes sein 🙂

    Jepp, mach ich.

    0
  • G

    headmyshoulder schrieb:

    GPC schrieb:

    Also ich habe gerade eben erfolgreich eine Verbindung auf den server hinbekommen. Vermutlich hast du dich vertippt 😉

    Nimmst Du die Daten von hier? Und wenn ja, normales FTP ohne verschlüsselung? Und welcher Port? Vielleicht ist mein Client auch nicht so intelligent (filezilla).

    Jappa, nehm ich. Und standardport. Hab's mit dem ftp-Konsolenprogramm gemacht 🤡

    gpc@desktop:~$ ftp c-plusplus.net
Connected to c-plusplus.net.
220 (vsFTPd 2.0.7)
Name (c-plusplus.net:gpc): magazin_public
331 Please specify the password.
Password:
230 Login successful.
Remote system type is UNIX.
Using binary mode to transfer files.
ftp> ls
200 PORT command successful. Consider using PASV.
150 Here comes the directory listing.
drwxr-sr-x    2 1003     33           4096 Dec 21  2007 Compiler
drwxrwsr-x    2 33       33           4096 Oct 03  2005 Debug_VC6
[.....]

    Kann's sein dass du ein O (Buchstabe) anstatt der Ziffer 0 nimmst (direkt nach pE1)?

    0
  • H

    GPC schrieb:

    Jappa, nehm ich. Und standardport. Hab's mit dem ftp-Konsolenprogramm gemacht 🤡

    Danke, mit dem Konsolenprogramm komm ich jetzt rein. Allerdings hab ich jetzt das Problem das ich kein Verzeichnis anlegen kann:

    mkdir template_metaprogrammierung_numerik

    liefert

    550 Create directory operation failed.

    Ansonsten seh ich alle Verzeichnisse und kann mir auch den Inhalt ansehen. Sry, falls ich mich zu doof anstelle.

    0
  • G

    headmyshoulder schrieb:

    GPC schrieb:

    Jappa, nehm ich. Und standardport. Hab's mit dem ftp-Konsolenprogramm gemacht 🤡

    Danke, mit dem Konsolenprogramm komm ich jetzt rein. Allerdings hab ich jetzt das Problem das ich kein Verzeichnis anlegen kann:

    mkdir template_metaprogrammierung_numerik

    liefert

    550 Create directory operation failed.

    Ansonsten seh ich alle Verzeichnisse und kann mir auch den Inhalt ansehen. Sry, falls ich mich zu doof anstelle.

    Sorry für die späte Antwort, die Woche war etwas stressig. Man kann atm tatsächlich kein Unterverzeichnis anlegen, habe da aber intern mal nachgefragt und sobald das wieder geht, geb ich dir hier Bescheid 😉

    0
  • G

    Okay, also die Berechtigungen sind jetzt frei. Ich hab das Verzeichnis schon mal unter dem Namen "template_metaprogrammierung_numerik" angelegt, du kannst jetzt die Inhalte selbst hochladen 😉

    0
  • G

    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 Compiliervorganges ausgewertet bzw. aufgefü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. Ausserdem 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 Templatemetaprgrammierung 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- und algorithmen und Compile-time Sequenzen. 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 add_const , remove_const zur Verfügung gestellt. Ein Teil dieser Funktionalität wird auch Einzug in den kommenden C++ Kompilerstandard finden.

    Zur Entwicklung von Meta-Algorithmen werden wir hier ausschliesslich Boost.MPL verwendet. Diese Bibliothek ist in ihrem Umfang und Anwendungsmöglichkeit 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 eine guter, wenn auch nicht einfacher Einstieg. Dort werden sogut 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 detailiert 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 Fliesskommazahlen 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. Ausserdem 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 die Operation (+,-,*,/) 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 Voila, man sieht dass mpl::_2::value nicht richtig ausgelö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

    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 grossen 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 klassischen 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] )

    Wie 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 ein alternative Implementation vorgestellt, welche sich durch eine einfache Konstrukion der Koeffizienten auszeichnet.

    Die Koeffizienten eines Verfahrens werden in mpl::vector 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 selber ist eine templatisierte Klasse, welcher 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ächsten überprüfen wir, ob butcher_a , butcher_b and 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 (ausser 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. Auch die letzt Stufe unterscheidet sich, 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 selber 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 Verfahren 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 Solver 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 Implementation

    Ein Nachteil der obigen Implementierung des Compile-Zeit 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 wird.

    Die 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 . Und calculate_stage übernimmt die Berechnung einer einzelnen Stage.

    Performance

    Die Performance der beiden Algorithmen ist etwa äquivalent. Ausserdem 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. Ausserdem stehen diese Algorithmen den ausgeschriebenen Algorithmen in Performance in nichts nach. Man kann diese Konzept natürlich auch einfach auf die Stepper mit Fehlerberechnung oder auf Stepper mit kontinuierlichen 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 weggemacht.

    Selbstverständlich erheben die beiden Version des Templatemetaprogrammierungs-Runge-Kutta-Verfahrens keinen Anspruch auf Allgemeingültigkeit. Sie wurde im Rahmen vom 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.

    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)

    0
Anmelden zum Antworten