Boost.uBLAS で行列式 - lapack版

boost C++ lapack uBLAS

今回は前回実装した行列式を lapack を用いて実装してみます．
前回と同様にまずLU分解をして，その結果を基に行列式を計算する流れになります．
LU分解は前回と同じく lapack.dgetrf を用います．その結果である上三角行列の対角成分の積を計算するだけで行列式を計算できます（ただし，行列式の性質から，置換が行われていた場合符号を変えて乗算します）．
以下に私の実装例を示します．

math.hpp

#ifndef MATH_HPP_20081015
#define MATH_HPP_20081015

#if defined(_MSC_VER) && (_MSC_VER >= 1020)
# pragma once
#endif


#include <boost/numeric/ublas/fwd.hpp>


namespace math {
    double determinant(const boost::numeric::ublas::matrix<double>& m);
}


#endif

math.cpp

#include "math.hpp"
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#define BOOST_NUMERIC_BINDINGS_USE_CLAPACK
#include <boost/numeric/bindings/lapack/gesv.hpp>
#include <boost/numeric/bindings/traits/ublas_matrix.hpp>
#include <boost/numeric/bindings/traits/ublas_vector.hpp>
#undef  BOOST_NUMERIC_BINDINGS_USE_CLAPACK


double math::determinant(const boost::numeric::ublas::matrix<double>& m)
{
    namespace ublas = boost::numeric::ublas;

    BOOST_UBLAS_CHECK(m.size1() == m.size2(), ublas::external_logic());

    ublas::permutation_matrix<int>             piv(m.size1());
    ublas::matrix<double, ublas::column_major> cm(m);
    int                                        info;

    info = boost::numeric::bindings::lapack::getrf(cm, piv);
    BOOST_UBLAS_CHECK(info >= 0, ublas::internal_logic());

    if (info == 0) {
        double det = 1;

        for (int i = 0; i < static_cast<int>(piv.size()); ++i) {
            det *= (i + 1 == piv(i)) ? +cm(i, i) : -cm(i, i);
        }

        return det;
    }
    else {
        return 0;
    }
}

math.cpp に一般行列の行列式を計算する関数 determinant を実装しました．math.cpp をプロジェクトに加え（他のソースコードと一緒にコンパイル），math.hpp を include することにより determinant が使用可能になります．
コードの簡単な説明ですが，bindings.lapack.getrf を呼び出すと，第一引数の下三角部分に単位下三角行列が，第一引数の上三角部分に上三角行列が，第二引数に置換行列が格納されます．その結果入力行列 A は PA = LU と分解されます．det(L) は明らかに 1 なので，det(A) は置換行列を参考に符号を変えながら U の対角成分の積を取っていけばよいことになります．lapack が返す置換行列のインデックスは uBLAS と違い 1 から始まっていることに注意してください（uBLAS でLU分解した場合，置換行列のインデックスは 0 から始まります）．
また，bindings.lapack.getrf が正の値を返した場合(info > 0 の場合)，LU分解後の上三角行列 U の戻り値に対応する対角成分が 0 になっていることを示します．このような場合行列式は明らかに 0 になるため，いちいち計算せずに 0 を返しています．


以下のテストコードで determinant の動作を確認できます．

test.cpp

#include "math.hpp"
#include <boost/numeric/ublas/matrix.hpp>


#ifdef _MSC_VER
# pragma comment(lib, "libf2c.lib")
# pragma comment(lib, "BLAS.lib")
# pragma comment(lib, "clapack.lib")
#endif


int main()
{
    boost::numeric::ublas::matrix<double> m(3, 3);

    m(0, 0) = -2; m(0, 1) = 2; m(0, 2) = -3;
    m(1, 0) = -1; m(1, 1) = 1; m(1, 2) = 3;
    m(2, 0) = 2;  m(2, 1) = 0; m(2, 2) = -1;

    // 計算結果を http://en.wikipedia.org/wiki/Determinant#Example と比較してみましょう
    std::cout << math::determinant(m) << std::endl;

    return 0;
}
  • VC9
  • g++ (GCC) 3.4.4 (cygming special)

でのコンパイルを確認しています．
実行してみると，正しく行列式を計算できていることが分かります．
エラーが出てコンパイルできない場合は，前回の記事を見て正しくインストールできているか確認してみてください．


次回はコレスキー分解を実装する予定です．
それでは．