メモの日々

■ [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";
    }
}