科学計算のためのc++テンプレートプログラミング

行列式の足し算などを行うとき、次のようなコードをかくと早いらしい。
処理によってはfortran, blasとかMKLより速いと公式のページにある。

ここのスライドを参考につくってみた。

 http://www.google.co.jp/url?sa=t&rct=j&q=advanced%20template%20programming%20&source=web&cd=1&ved=0CDUQFjAA&url=http%3A%2F%2Fitee.uq.edu.au%2F~conrad%2Fmisc%2Fsanderson_templates_lecture_uqcomp7305.pdf&ei=Xx4rUarGLvHsmAXbaw&usg=AFQjCNFiGSmasUA2pfRuHGdVExcyw5tUaA&bvm=bv.42768644,d.dGY&cad=rja

なぜ早いかというと足し算のあとに一時的な行列の実態が作られないからだそうだ。　バイナリーになおすとoperator+( , )から呼び出されているはずのGlueのコンストラクタがないらしい。

 こんな感じのmeta programmingした数学ライブラリーがarmadilloとして配布されている

#include<iostream>

class Glue;

class Matrix{
public:
    Matrix();
    Matrix( int in_rows, int in_cols );
    void set_size( int in_rows, int in_cols );
    void clear();
    void printm();

    Matrix( const Matrix & X );
    const Matrix & operator=( const Matrix & X );

    Matrix( const Glue & X );
    const Matrix & operator=( const Glue & X );

    int rows;
    int cols;
    double *data;
};

class Glue{
public:
    const Matrix& A;
    const Matrix& B;
    Glue( const Matrix& in_A, const Matrix & in_B )
    : A( in_A )
    , B( in_B ){}
};

/*--------------------------
    Glue
---------------------------*/

const Glue operator+( const Matrix& A, const Matrix& B ){
    return Glue(A,B);   
}

/*--------------------------
    Matrix
---------------------------*/
Matrix::Matrix(){
}

Matrix::Matrix( int in_rows, int in_cols){
    cols = in_cols;
    rows = in_rows;
    set_size( rows, cols );
}

void Matrix::set_size( int in_rows, int in_cols ){
    data = new double [ in_rows*in_cols ];   
    for( int i=0; i<in_rows*in_cols; i++ ) data[i] = i;
}

void Matrix::clear(){
    for( int i=0; i<rows*cols; i++ ) data[i] = 0.0;
}

void Matrix::printm( ){
    for( int j=0; j<rows; j++ ){
        for( int i=0; i<cols; i++ ){
            std::cout << " " << data[i+j*cols];
        }
        std::cout  << std::endl;
    }

}

const Matrix & Matrix::operator=( const Glue & X ){
    const Matrix & A = X.A;
    const Matrix & B = X.B;

    set_size( A.rows, A.cols );

    for( int i=0; i < A.rows * A.cols; i++ ){
        data[i] = A.data[i] + B.data[i];
    }
    return *this;
}

Matrix::Matrix( const Glue & X ){
    operator=(X);
}

int main(){
    Matrix m1(4,4);
    Matrix m2(4,4);
    Matrix m3(4,4);

    m3.clear();

    m3 = m1 + m2;

    //m1.printm();
    //m2.printm();
    //m3.printm();
return 0;}

mainのなかのをコメントアウトすれば結果が出てくるよ。 実際OpenFoamとかはこういう書き方になっているが、Openfoamのは一時ファイルを作るやり方で書いてある。

OpenFOAMをmetaできるか考えてみたが、privateの部分をpublicに治す必要がある。　つまりクラスを使う意味がなくなってしまうのでは？ そして、やっぱりコードの抽象度は実行速度とは反比例する とおもう。情報の人はこのように書きたいのだろうけど、流体コード屋としては非常に書きづらい。　もっと低いレベルで書いたほうがしっくりくる。