Boost.uBLAS で逆行列 - lapack版
今回は前回実装した逆行列を lapack を用いて実装してみます.
前回と同様にまずLU分解をして,その結果を基に逆行列を計算する流れになります.
LU分解は前回と同じく lapack.dgetrf を用います.その結果を lapack.dgetri に渡すだけで逆行列を計算できます.
LU分解のときのように行列を一生懸命切り出す必要がないため,非常に簡単に実装できうれしい限りです.
以下に私の実装例を示します.
math.hpp
#ifndef MATH_HPP_20081014 #define MATH_HPP_20081014 #if defined(_MSC_VER) && (_MSC_VER >= 1020) # pragma once #endif #include <boost/numeric/ublas/fwd.hpp> namespace math { void invert(const boost::numeric::ublas::matrix<double>& m, boost::numeric::ublas::matrix<double>& mi); } #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 void math::invert(const boost::numeric::ublas::matrix<double>& m, boost::numeric::ublas::matrix<double>& mi) { 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()); info = boost::numeric::bindings::lapack::getri(cm, piv); BOOST_UBLAS_CHECK(info == 0, ublas::internal_logic()); #if BOOST_UBLAS_TYPE_CHECK BOOST_UBLAS_CHECK( ublas::detail::expression_type_check( ublas::prod(m, cm), ublas::identity_matrix<double>(m.size1()) ), ublas::internal_logic() ); #endif mi = cm; }
math.cpp に一般行列の逆行列を計算する関数 invert を実装しました.デバッグ時には実際に計算結果を用いて m x mi が単位行列になるかを確認するコードを追加しています.math.cpp をプロジェクトに加え(他のソースコードと一緒にコンパイル),math.hpp を include することにより invert が使用可能になります.
コードの簡単な説明ですが,bindings.lapack.getrf を呼び出すと,第一引数の下三角部分に単位下三角行列が,第一引数の上三角部分に上三角行列が,第二引数に置換行列が格納されます.
しかしそんなことを気にする必要はなく,得られた二つの行列をそのまま単に bindings.lapack.getri に渡すことにより逆行列を計算することができます.
その結果を関数の引数 mi に保存して終了です.
以下のテストコードで invert の動作を確認できます.
test.cpp
#include "math.hpp" #include <boost/numeric/ublas/io.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() { namespace ublas = boost::numeric::ublas; ublas::matrix<double> m(3, 3), mi; m(0, 0) = 3; m(0, 1) = 2; m(0, 2) = 1; m(1, 0) = 1; m(1, 1) = 3; m(1, 2) = 1; m(2, 0) = 2; m(2, 1) = 2; m(2, 2) = 1; math::invert(m, mi); // 以下結果の確認 std::cout << mi << std::endl; BOOST_UBLAS_CHECK(ublas::detail::expression_type_check(ublas::prod(m, mi), ublas::identity_matrix<double>(3)), ublas::internal_logic()); BOOST_UBLAS_CHECK(ublas::detail::expression_type_check(ublas::prod(mi, m), ublas::identity_matrix<double>(3)), ublas::internal_logic()); return 0; }
- VC9
- g++ (GCC) 3.4.4 (cygming special)
でのコンパイルを確認しています.
実行してみると,正しく逆行列を計算できていることが分かります.
エラーが出てコンパイルできない場合は,前回の記事を見て正しくインストールできているか確認してみてください.