B級科学者もどきの憂鬱

とある理系になりきれない奴のつれづれなる活動記

スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

行列計算ライブラリEigenで固有値の計算

お久しぶりです。もうすぐ年明けですね。
もう何ヶ月も全然書いてなかったので、
次の更新はTIPSの話にしたかったのですが、
記事にできるほど進んでません!
というわけで今回は全然違うテーマです。

最近は数値シミュレーションを色々とやっています。
その過程で、行列計算にEigenというC++用ライブラリを導入しました。
色々見てる限り、結構速そうで導入も楽で使いやすそうです。

行列計算だけでなく、普通の配列を扱う計算で使っても、
アセンブラレベルで並列化したりしてくれるので、
普通にC++で書くより速いそうです。
これに関してはまだ確かめてませんけど。

導入は、公式サイトから圧縮ファイルをダウンロードし、
解凍して出来たファイルをインクルードパスに追加するだけです。
ヘッダファイルだけで出来ているので、コンパイル等は不要です。

例として、実数行列の固有値を求めるコードを書いてみます。

#include <stdio.h>
#include <Eigen/Dense>  // 固有値計算ルーチンの入ったインクルードファイル

int main() {
  Eigen::Matrix<double, 3, 3> A;  // 実数行列
  A << 1,2,3,4,5,6,7,8,9;  // Aの行列要素を代入

  Eigen::EigenSolver< Eigen::Matrix<double, 3, 3> > s(A);
  std::cout << "固有値\n" << s.eigenvalues() << std::endl;
  std::cout << "固有ベクトル\n" << s.eigenvectors() << std::endl;

  return 0;
}

EigenSolverが固有値を求めるテンプレートクラスで、
テンプレート引数に、求めたい行列の型名を入れて使います。
行列型の型名が割と長くなるので、typedef等で短い名前に置き換えた方が、
間違いが少なくなると思います。

複素行列の固有値を求めたい場合は、Aの型名の、
Eigen::Matrix<double, 3, 3>
doublestd::complex<double>に置き換えて、
(もちろん、EigenSolverのテンプレート引数も一緒に置き換えて)
Eigen::Matrix<std::complex<double>, 3, 3>
とすれば良い、だけではありません

これだけだと、EigenSolverのコンパイルに失敗します。
少なくとも、私が使った環境ではエラーになります。
(Eigen3.2.0、Visual Studio Express 2012 for Windows Desktop)
私は最初、これにハマって数十分悩みました。
というかこれが、今回の記事を書いた動機です。

解決法は、EigenSolverの代わりに、ComplexEigenSolverを使うことです。
これでエラーは出なくなります。

さらに詳しく知りたい方は、公式リファレンスを参照して下さい。
スポンサーサイト

#defineでシングルトン化

デザインパターンの中で、おそらく最もシンプルなのが、
シングルトンパターンでしょう。
他のパターンはあんまり知らないけど、
これだけは使えるという人も多いです。

このパターン、C++のプログラミング系のサイトで、
よくテンプレート化出来ないか取り上げられてますよね。

ただ、シングルトンパターンをテンプレートの継承で行う場合、
どうやってもコンストラクタがprotectedになり、
完全なシングルトンにすることが出来ません。
継承を使わなくとも、色々と問題点が多く、
シンプルな解は中々無いようです。

じゃあテンプレートじゃなくて、マクロでやればいいんじゃね?
と思って作ってみました。

#define SINGLETONIZE(ClassName) \
    private: \
        ClassName(); \
        ~ClassName(); \
        ClassName(const ClassName&); \
        ClassName& operator=(const ClassName&); \
    public: \
        static ClassName& GetInstance() { \
            static ClassName object; \
            return object; \
        }

とても簡単に出来ちゃいました。

シングルトンパターンを知っている人なら、
ほとんど説明はいらないと思いますが、
一応、使い方はこんな感じになります。
シングルトン化したいクラスを、Singleとします。

class Single {

    // メンバ色々……

    SINGLETONIZE(Single)

    // メンバ色々……

};

という風に、クラス宣言内に書けばOKです。

シングルトン化したいクラス名以外のものをマクロの引数に渡すと、
多分全ての場合でコンパイラがエラーを吐きますので、
マクロでよくある様々な障害とも無縁だと思います。

スレッドセーフな実装ではありませんので、
マルチスレッドでは直接使えませんが、
多少改変すれば出来るはずです。

こんなマクロを作っている人を案外見かけないのですが、
どうしてなのでしょうか?
あまりにも自明な解過ぎて載せる価値が無いという事か、
はたまた、マクロなんて邪道ってことなんでしょうかね。
何事も適材適所だと思いますが。

基底クラスの後付け

ちょっと開発中に思いついたことの覚書も兼ねて、
久々にプログラミングカテゴリの記事です。

C++において、ある二つのクラスに、
同じ働きをするメソッドがありますが、
それらは、そのメソッドを実装した、
同じ基底クラスを継承していないとします。

このとき、既存のクラスに変更を加えること無く、
仮想的に同じクラスを継承しているように
見えるようにするにはどうすれば良いか、
というのが今回の話です。

上の説明だけでは抽象的で分かりにくいと
思いますので、具体的にコードを書いてみます。

struct Class1 {
    void Run();
};

struct Class2 {
    void Run();
};

見ての通り、同じ名称のコードが宣言されています。
処理内容もほぼ同じだとします。

これらには、基底クラスが無いので、
例えばこんなことは出来ません。

struct Base {
    virtual void Run()=0;
    virtual ~Base() {}
};

Base *object[2];
object[0]=new Class1();
object[1]=new Class2();
...

これを何とか出来るようにするには
どうすれば良いか、という話です。

前置きが長くなりましたが、解決方法は割と単純です。
次のようなテンプレートクラスを用意します。

template < class T > class Adapter : public Base {
private:
    T object;
public:
    virtual void Run() {
        object.Run();
    }
};

これを使って、先程の例を次のように変更します。

Base *object[2];
object[0]=new Adapter<Class1>();
object[1]=new Adapter<Class2>();
...

これで、既存のクラスに一切変更を加えずに、
基底クラスのようなものを作れました。

ところでこれ、テンプレートが使えない場合、
どうすればいいんでしょうね。
さっぱり見当が付きません。


>>管理者用コメントの方
ありがとうございます。
最近サボるということを覚えてきたので、
あまり無理をしない程度に頑張ります。

DirectXは何故floatばかり使う?

もはや一月更新が板に付いてきました。
書く時間自体は多少あったのですが、
如何せん中々文章がまとまりませんでした。
結局3つぐらい下書き状態で放置です。
今回は、DirectXの話です。

DirectXのプログラムを書いていると、
座標などがほとんどfloatで扱われていることに気付きます。
何故、doubleを使わないでfloatなのでしょう?
CPUは大抵、doubleの方が計算速いのに……。
というわけで、ちょっと調べてみました。


まず一つ目の理由。それは、計算が速いからです。
あれ、さっきdoubleの方が計算速いとか言ってなかったっけ?
いえいえ、それはあくまでCPUの話です。

DirectXはグラフィックライブラリなので、
ほとんどの計算処理を行うのはGPUです。
そして、GPUは大抵floatの方が計算が速いのです。


二つ目の理由は、メモリ使用量が少ないからです。
まぁ、doubleに比べればサイズは半分ですからね。
ですが、単純にそれだけじゃありません。

DirectX7までは、3D頂点の情報を入れるための専用の構造体がありました。
これには普段使わないメンバも多かったので、メモリの無駄を省く為、
DirectX8からは、プログラマが必要なメンバだけを
定義した構造体が使えるようになりました。

DirectX側には、どんな風にメンバを宣言したかを伝えておけばOKです。
例えば、最初の12byteは、float型でx,y,z座標が入っている、
次の4byteに色の情報が入っている、などです。

ここで、頂点の座標等にfloatではなくdoubleを使うとどうなるでしょうか。
もし、doubleが頂点情報を格納する構造体にあると、
4byteのメンバと8byteのメンバが混在することになるので、
コンパイラにより、パディングが入ってしまう可能性があります。

そうなると、本来のメンバに必要な領域だけでなく、
一切使われない領域が出来てしまう、ということが有り得るのです。
メモリは出来る限り節約したいので、これは避けたいですよね。

つまり、頂点構造体のメンバを4byteに統一するために、
doubleではなくfloatを使っているわけです。
これならパディングは一切入りません。


二つ理由を紹介しました。
実際には、これらは相互に関連しているのでしょう。
floatがよく使われるから、GPUの計算もfloatの方を速くし、
計算が速いから、ますますfloatばかり使われる……と。
まぁこれは私の推測でしか無いですが。

ちなみに、最近はGPUで数値計算を行わせることも多いので、
double型などの計算もかなり速くなっているようです。
それでも、floatの方が数倍速いのですが。

多重ループを抜けるための6つの方法

結局、6月中には更新できなかったどころか、
また一ヶ月ほど空いてしまいました。

最近、空いた時間というのは、
いつの間にかあるものじゃなくて、自分で作るものである、
ということを半ば強制的に実感させられています。
7月過ぎたら多少は時間が作れると信じてる!

今回は久々にプログラミングの話です。
まぁ、備忘録的なものですが。

二重三重にも積み重なったループを、
一気に抜けたい時ってありますよね。
ですが、私が知っている限り、プログラミング言語で、
多重ループを抜けるための言語仕様があるものはありません。

これを実現するためには様々な方法がありますが、
代表的と思われるものを6つまとめてみました。
言語は全てC++ですが、他の言語でも応用できるでしょう。

以下のソースコードは、xMax, yMax, zMaxというものが
どこかで宣言されていると思って下さい。
もしかしたらglobal変数かも知れませんし、
#defineで定義された定数かも知れません。


1.gotoを使う

int x, y, z;
for(x=0; x<xMax; x++) {
    for(y=0; y<yMax; y++) {
        for(z=0; z<zMax; z++) {
            if(終了条件) goto for_exit;
            何かしらの処理
        }
    }
}
for_exit:

一番シンプルな方法です。
ループがいくら重なっても同じやり方で出来ますし、
処理コストも非常に少ないです。
たまにgotoを使うのを極端に嫌う人がいますが、
何事も適材適所だと思います。


2.脱出用フラグを使う

int x, y, z;
bool exitFlag=false;
for(x=0; x<xMax; x++) {
    for(y=0; y<yMax; y++) {
        for(z=0; z<zMax; z++) {
            if(終了条件) {
                exitFlag=true;
                break;
            }
            何かしらの処理
        }
        if(exitFlag) break;
    }
    if(exitFlag) break;
}

割とメジャーな方法だと思います。
重ねたループの数だけフラグの判定をする必要があるので、
見た目は多少煩雑になり、処理コストも少し掛かります。

ですが、一番内側のループ以外にも処理がある場合、
例えば上記の場合、zのループ終了後にyループ内で別の処理がある場合、
フラグ判定文の位置を変えるだけで対応できます。
gotoで同じ事をやろうと思うと、スパゲッティになります。


3.ループ終了条件を無理矢理満たすようにする

int x, y, z;
for(x=0; x<xMax; x++) {
    for(y=0; y<yMax; y++) {
        for(z=0; z<zMax; z++) {
            if(終了条件) {
                x=xMax-1;
                y=yMax-1;
                z=zMax-1;
            }
            何かしらの処理
        }
    }
}

脱出用フラグを使う代わりに、
ループ変数自体を変更してしまうものです。
これなら脱出用フラグを使う必要はありません。

ですが、脱出のためにいくつも代入文を書く羽目になるので、
途中で間違える確率は上がります。
ループ上限回数を変更したい場合も、
複数箇所変更しなければならない可能性があります。

ただ、内側のループ2つだけ抜けるのと、3つ全部抜けるのを、
条件によって切り替えたい、などの場合は、
こちらの方法が使えます。

脱出用フラグで同じ事をしようとすると、
抜け方によってそれぞれフラグを用意しなければならず、
ミスも増える可能性が高いです。

それと、if文内部で、ループ変数を上限値から1引いた値にしていますが、
引かなくても大抵の場合は動作します。

しかし、例えばxMaxがINT_MAXだった場合、
x=xMax;など、1を引かない形にしてしまうと、
次のxのループが始まるとき、再設定式x++が行われるので、
桁あふれを起こして、xが負の値になります。

そうなると、継続条件式がtrueになるため、
意図通りにループを脱出できません。


4.ループ部分だけ関数にする

void Loop(void) {
    int x, y, z;
    for(x=0; x<xMax; x++) {
        for(y=0; y<yMax; y++) {
            for(z=0; z<zMax; z++) {
                if(終了条件) return;
                何かしらの処理
            }
        }
    }
}

ループを丸ごと関数にして、
returnで一気に抜ける、という方法です。
gotoを使った場合と同じような特徴があるにも関わらず、
gotoを使わないで済む上に、構造化も図れます。

ただ、このようなループがプログラム内に多数存在する場合、
関数の名前のストックが無くなってくるという、
ものすごく現実的な問題があります。
まぁ、私が直面しているだけかも知れませんが。

適当に名付けるとプログラムの流れが見えにくくなるので、
出来るだけ、何をやる関数なのか分かる名前で作りましょう。


5.例外を利用する

try {
    int x, y, z;
    for(x=0; x<xMax; x++) {
        for(y=0; y<yMax; y++) {
            for(z=0; z<zMax; z++) {
                if(終了条件) throw(-1);
                何かしらの処理
            }
        }
    }
} catch (int) {
    ループ終了時の処理
}

例外を投げた直後に、catchブロックまで
制御が飛んでいくという性質を利用したものです。
これも1.や4.とほぼ同じ特徴を持ちます。

例外をループ脱出に使うのに抵抗のある人もいるかもしれませんが、
プログラムというのは、誰にでも見やすく書けるのなら
基本的に何でもありです。

投げる例外は別に何でもいいですが、
どこかでExitLoop型等を宣言しておいて、
それを投げることにすると、例外の目的がはっきりします。


6.無理矢理ループを一重にする

int x, y, z, xyz, xyzMax;
xyzMax=xMax*yMax*zMax;
for(xyz=0; xyz<xyzMax; xyz++) {
    x=xyz/yMax/zMax;
    y=xyz/zMax%yMax;
    z=xyz%zMax;
    if(終了条件) break;
    何かしらの処理
}

無理矢理一つのループにして、そのループ変数の値から、
本来のループ変数の値を復元する、という力技な方法です。
今までの例の中で一番柔軟性が無いと言えるかも知れません。

四重以上のループとなると、
ループ変数の復元が複雑になるのでまず行われませんが、
二重ループ程度ではそんなに手間ではありません。

こんな強引なことを行うのは、二次元配列として扱うべきものを
プログラム上は一次元配列で格納している、等の特殊な場合が多いです。
行列計算ライブラリなどでたまに見かける実装です。

この方法には、たとえxMax, yMax, zMaxがINT_MAX以下でも、
xyzMaxがINT_MAXを超える可能性があるという問題もあります。
なので、あまりループ回数の多いものは作れません。


以上、6つの方法を紹介しました。
これらの派生形も色々と考えられますので、
形にとらわれず、目的に応じて工夫して使って下さい。
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。