Boost.uBLAS で逆行列 - lapack版

boost C++ lapack uBLAS

今回は前回実装した逆行列を 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)

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


次回は lapack を用いた行列式を実装する予定です．
それでは．