Boost.uBLAS でLU分解

boost C++ uBLAS

今回はLU分解を実装します．LU分解を計算する関数はもちろん lapack にもあるのですが，それは次回紹介するとして，今回は uBLAS に実装されている lu_factorize を用いてLU分解してみようと思います．uBLAS の lu_factorize の結果を三角行列に切り出すだけなので非常に簡単です．以下に私の実装例を示します．

math.hpp

#ifndef MATH_HPP_20081012
#define MATH_HPP_20081012

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


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


namespace math {
    void lu_factorize(boost::numeric::ublas::matrix<double> A, boost::numeric::ublas::triangular_matrix<double, boost::numeric::ublas::unit_lower>& L, boost::numeric::ublas::triangular_matrix<double, boost::numeric::ublas::upper>& U);
}


#endif

math.cpp

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


void math::lu_factorize(boost::numeric::ublas::matrix<double> A, boost::numeric::ublas::triangular_matrix<double, boost::numeric::ublas::unit_lower>& L, boost::numeric::ublas::triangular_matrix<double, boost::numeric::ublas::upper>& U)
{
    namespace ublas = boost::numeric::ublas;

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

    ublas::lu_factorize(A);

    ublas::triangular_matrix<double, ublas::unit_lower> CL = ublas::triangular_adaptor<ublas::matrix<double>, ublas::unit_lower>(A);
    ublas::triangular_matrix<double, ublas::upper> CU = ublas::triangular_adaptor<ublas::matrix<double>, ublas::upper>(A);

    L.assign_temporary(CL);
    U.assign_temporary(CU);
}

math.cpp にLU分解を計算する関数 lu_factorize を実装しました．
関数 lu_factorize は M 次正方行列 A を受け取り，これを M 次単位下三角行列 L と M 次上三角行列 U の積に因数分解します．M 次単位下三角行列が関数の引数 L に，M 次上三角行列が関数の引数 U に保存されます．
lu_factorize の第一引数が const 参照になっていないことに注意してください．こうすることでコンパイラはより効率的なコードを生成することができるかもしれません（しかし Meyers が指摘しているように，このコードは賢さを求めすぎ，明白さを失ってしまったかもしれません）．しかし，BOOST_UBLAS_CHECK の前にコピーが発生するのを許容できない方は const 参照で受け取りましょう（私がそうです．今回は非参照で受け取れば効率のよいコードを生成できる可能性があるということを言いたかっただけです）．
math.cpp をプロジェクトに加え（他のソースコードと一緒にコンパイル），math.hpp を include することにより lu_factorize が使用可能になります．


以下のテストコードで lu_factorize の動作を確認してみます．

test.cpp

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


int main()
{
    namespace ublas = boost::numeric::ublas;

    ublas::matrix<double> A(2, 2);

    A(0, 0) = 4; A(0, 1) = 3;
    A(1, 0) = 6; A(1, 1) = 3;

    ublas::triangular_matrix<double, ublas::unit_lower> L;
    ublas::triangular_matrix<double, ublas::upper>      U;

    math::lu_factorize(A, L, U);

    // 以下結果の表示
    // 計算結果を http://en.wikipedia.org/wiki/LU_decomposition#Small_Example と比較してみましょう
    std::cout << L << std::endl
              << U << std::endl
              << ublas::prod(L, U) << std::endl;

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

でのコンパイルを確認しています．
実行してみると，正しくLU分解できていることが分かります．


次回は lapack を用いたLU分解を実装する予定です．
それでは．