メモの日々

■ [c][math] LAPACKを使って行列の固有値計算

線形代数計算ライブラリのLAPACKを使ったのでちょっとだけメモ。

LAPACKを説明したサイトは色々あるけど、数値計算マニュアル：線型代数 (FAQ)がよくまとまっていてよかった。日本計算工学会誌「計算工学」からリンクされている「BLAS, LAPACKチュートリアル パート1」（Vol16,No.2）、「BLAS, LAPACKチュートリアル パート2」（Vol.16,No.3）も参考になる。

実対称行列の固有値と固有ベクトルを計算してくれる関数 dsyev_() を使ったサンプルを書いてみた。関数の戻り値が何を表すのかが不明。

#include <cmath>
#include <iostream>
#include <stdexcept>
#include <vector>

using std::vector;

// f2c.h, clapack.h から必要な宣言を抜粋
extern "C" {
    typedef long int integer;
    typedef double doublereal;

    int dsyev_(
            char *jobz,
            char *uplo,
            integer *n,
            doublereal *a,
            integer *lda,
            doublereal *w,
            doublereal *work,
            integer *lwork,
            integer *info);
}

/**
 * 実対称行列の固有値/固有ベクトルを計算する
 *
 * @param matrix 計算対象の行列。上三角部分だけが設定されていればよい。
 *               matrixの内容は破壊される(固有ベクトルが設定される)。
 *
 * @return matrixの全固有値からなるvectorと、全固有ベクトルのvector。
 *         戻り値の0番目のvectorが全固有値、
 *         戻り値の1番目のvectorが0番目の固有値に対する固有ベクトル、
 *         戻り値の2番目のvectorが1番目の固有値に対する固有ベクトル、
 *         ...
 */
std::vector<std::vector<double> > eigen(std::vector<double>& matrix) {
    if (matrix.empty()) return vector<vector<double> >();

    char jobz = 'V';                // 固有ベクトルも計算する
    char uplo ='U';                 // 上三角部分を使用。ただし列優先。
    integer n = static_cast<integer>(std::sqrt(matrix.size())); // matrixの行数
    integer lda = n;                // matrixのleading dimension
    vector<doublereal> work(n * 3); // dsyev_()の作業領域
    integer lwork = static_cast<integer>(work.size());          // workのサイズ
    integer info = 0;               // dsyev_()の処理結果

    vector<double> eigenvalues(n);
    dsyev_( &jobz,
            &uplo,
            &n,
            &matrix[0],
            &lda,
            &eigenvalues[0],
            &work[0],
            &lwork,
            &info);

    if (info) {
        throw std::runtime_error("error: dsyev_()");
    }

    vector<vector<double> > result(n + 1, vector<double>(n));
    result[0] = eigenvalues;

    // LAPACKは列優先
    for (integer col = 0; col < n; ++col) {
        for (integer row = 0; row < n; ++row) {
            result[col + 1][row] = matrix[col * n + row];
        }
    }
    return result;
}

int main() {
    vector<double> m;
    m.push_back(0);
    m.push_back(1); // ここの値は使われない
    m.push_back(1);
    m.push_back(-1);

    vector<vector<double> > eigens = eigen(m);

    if (eigens.empty()) {
        throw std::runtime_error("BUG: no eivenvalue.");
    }

    std::cout << "eigenvalues:";
    for (std::size_t i = 0; i < eigens[0].size(); ++i) {
        std::cout << " " << eigens[0][i];
    }
    std::cout << "\n";

    for (std::size_t i = 1; i < eigens.size(); ++i) {
        std::cout << "eigenvector" << i << ":";
        for (std::size_t j = 0; j < eigens[i].size(); ++j) {
            std::cout << " " << eigens[i][j];
        }
        std::cout << "\n";
    }
}

■ [life] 健康診断を受けた2011

今年も受診できた。去年より混んでいる印象。老人だらけで、そりゃ国保の保険料高くなるわ…。

まず身長体重。柱に寄り掛かっては駄目と言われた。機械が変わったのかもしれない。次血圧。そしてしばらく待つ。いつものパターン。

待っている間に、あるおばあさんの名前が呼ばれ、それに対し間違ったおばあさんが説明を聞いて次の検査へ行ってしまうというトラブルが発生した。レントゲンの方で名前が一致しないことに気付き解決した模様。

おれの番も来て、レントゲン。名前を呼ばれる前にレントゲン室の扉を開けてしまったのは失敗だった気がする。でも上着を脱いだら中へ入れと言われたんだよなあ。続けて尿検査、採血、心電図。しばらく待って医師による問診。胴回りは78cm。去年と同じなのか。

今年は1時間半掛かった。料金は500円。後日便を採取して提出しないといけないのだ。未経験なので怖い。

■ [life] 歯医者へ行く 2回目

結局、今回も検査をされただけで治療行為はなし。

歯周ポケットの深さを、全歯について調べられた。所によりわりと痛い。衛生士さんがマイクで機械に「次」「戻れ」「3」などど指示しながら記録を録っていた。奥歯の方に深い歯周ポケットがある。

最後にレントゲン。メリーゴーランドのような音楽と共に機械が顔の回りをぐるぐる回っていた。

院長先生は不在なようで、若い先生が少しだけ顔を出してくれた。料金は3600円。

■ [life] 歯医者へ行く 3回目

前回撮ったレントゲンの写真を見る。下の親知らずが横向きに生えていた。が、痛みは無いのでそこはスルー。

続いて歯磨きについて色々質問される。答えた後歯の磨き方を教わったが、目新しい情報はない。

次に汚れが紫になる薬を塗られた自分の歯を鏡で見る。歯間に汚れが残っているのが分かる。歯ブラシを渡され、教わった磨き方で磨かされた。途中から衛生士さんに代わってもらった。

本日は以上で1700円。未だ治療に入らない。歯周病はたいしたことないのだな。

■ [shell][howto] 標準出力と標準エラーを入れ替える

• 36.16. n>&m: Swap Standard Output and Standard Error （Unix Power Tools, 3rd Edition）

に、

$ command 3>&2 2>&1 1>&3

というイディオムが紹介されていた。これでcommandの標準出力と標準エラーが入れ替わる。

動作確認。

% bash

$ ls
tmp/

$ ls tmp tmp2
ls: tmp2: そのようなファイルやディレクトリはありません
tmp:

$ ls tmp tmp2 3>&2 2>&1 1>&3 | sed 's/tmp/XXX/'
tmp:
ls: XXX2: そのようなファイルやディレクトリはありません

標準エラーに対してだけ置換が行われたのでOK。ちなみに、zsh 4.2.0（古い）だとうまくいかなかった。

■ [shell][howto] 標準エラーだけを操作する

上のようにすれば標準エラーだけをパイプに渡せるけど、標準出力と標準エラーが入れ替わってしまう。標準出力は維持し、標準エラーに対してだけ置換を行った結果を標準エラーに出力して欲しい。

上のようにした後でもう一度入れ替えればいいはず。

$ (ls tmp tmp2 3>&2 2>&1 1>&3 | sed 's/tmp/ttt/') 3>&2 2>&1 1>&3 | sed 's/t/X/'
ls: ttt2: そのようなファイルやディレクトリはありません
Xmp:

もっと簡単にやる方法あるのかな？

■ [life] 歯医者へ行く 4回目

今日は歯石を取ってもらった。やっと治療的なことをしてもらえたぞ。

口用の穴が開いたタオルを顔に掛けられて実施。冷却水（？）が結構飛び散ってきた。機械でカリカリやられたのだけれど、痛みはそれ程ない。たまに鈍痛があって、そのときは緊張した。薬を塗られて終了。歯が白くなっていた。2000円。