Boost.uBLAS で行列式
前回は逆行列を紹介しました.逆行列の次にくるものと言えば,やはり行列式でしょう(どちらかというと逆かもしれませんが).
逆行列のときと同様に,行列式を求める関数もuBLASには存在しないため,自分で作るなりどこかからもってくるなりしなくてはなりません.
しかし,ありがたいことに,行列式のコードも既に多くの人が書いてくれています.
- Ublas mailing page: Re: ublas determinant and norm
- Re: matrix expressions (was find2 required by Matrix)
実装は逆行列のときとほとんど同じです.
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
また,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 のインストールについて紹介しようと思います.
それでは.