Boost.uBLAS で行列式

前回は逆行列を紹介しました.逆行列の次にくるものと言えば,やはり行列式でしょう(どちらかというと逆かもしれませんが).
逆行列のときと同様に,行列式を求める関数もuBLASには存在しないため,自分で作るなりどこかからもってくるなりしなくてはなりません.
しかし,ありがたいことに,行列式のコードも既に多くの人が書いてくれています.

実装は逆行列のときとほとんど同じです.

math.hpp

#ifndef MATH_HPP_20080916
#define MATH_HPP_20080916

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


#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/matrix.hpp>


namespace math {
    template <class M> double determinant(const M& m);
}


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

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

    ublas::matrix<double>       lu(m);
    ublas::permutation_matrix<> pm(m.size1());

    ublas::lu_factorize(lu, pm);

    double det(1);

    typedef ublas::permutation_matrix<>::size_type size_type;

    for (size_type i = 0; i < pm.size(); ++i) {
        det *= (i == pm(i)) ? +lu(i, i) : -lu(i, i);
    }

    return det;
}


#endif

math.hpp に行列式を求める関数 determinant を実装しました.
math.hpp を include するだけで determinant が使用可能になります.
determinant の戻り値には少し注意が必要です.
もし戻り値を M::value_type にしてしまうと,matrix などで桁があふれてしまうおそれがあるためです(matrix行列式が short int に収まる保証はありません).
また,lu の要素が整数型ではLU分解を正確に計算できないため,結局内部では float や double あるいは long double で
計算することになります(そのため lu の型は M ではなく matrix になっています).
以上の理由により,今回の例では戻り値を double にしました.

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

test.cpp

#include "math.hpp"


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

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

    std::cout << math::determinant(m) << std::endl; // 18 になるはず

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

でのコンパイルを確認しています.
実行してみると,正しく行列式が求められていることが分かります.


次回は Boost.uBLAS+Boost.Bindings+lapackの利用に向け,
Bindings と lapack のインストールについて紹介しようと思います.
それでは.