Boost.uBLAS で一般係数行列の線形方程式系(System of linear equations)の解の計算

今回は AX = B の解行列 X あるいは Ax = b の解ベクトル x を計算する関数を実装します.uBLAS にはLU分解を計算する機能が既にあるため,これらの機能を簡単に実装することができます.
また,uBLAS による連立方程式の解き方は

で既に詳細に説明されており,私の能力ではこれを超える記事を書くことはできないのですが,これらの結果をまとめて誰でも使えるような関数に実装することには多少の意義があると思い,以下に私が実装した例を示します.

math.hpp

#ifndef MATH_HPP_20081105
#define MATH_HPP_20081105

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


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


namespace math {
    void solve(const boost::numeric::ublas::matrix<double>& A, const boost::numeric::ublas::vector<double>& b, boost::numeric::ublas::vector<double>& x);
    void solve(const boost::numeric::ublas::matrix<double>& A, const boost::numeric::ublas::matrix<double>& B, boost::numeric::ublas::matrix<double>& X);
}


#endif

math.cpp

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


void math::solve(const boost::numeric::ublas::matrix<double>& A, const boost::numeric::ublas::vector<double>& b, boost::numeric::ublas::vector<double>& x)
{
    namespace ublas = boost::numeric::ublas;

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

    ublas::matrix<double>       lhs(A);
    ublas::vector<double>       rhs(b);
    ublas::permutation_matrix<> pm(A.size1());

    ublas::lu_factorize(lhs, pm);
    ublas::lu_substitute(lhs, pm, rhs);

#if BOOST_UBLAS_TYPE_CHECK
    BOOST_UBLAS_CHECK(ublas::detail::expression_type_check(ublas::prod(A, rhs), b), ublas::internal_logic());
#endif

    x.assign_temporary(rhs);
}


void math::solve(const boost::numeric::ublas::matrix<double>& A, const boost::numeric::ublas::matrix<double>& B, boost::numeric::ublas::matrix<double>& X)
{
    namespace ublas = boost::numeric::ublas;

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

    ublas::matrix<double>       lhs(A);
    ublas::matrix<double>       rhs(B);
    ublas::permutation_matrix<> pm(A.size1());

    ublas::lu_factorize(lhs, pm);
    ublas::lu_substitute(lhs, pm, rhs);

#if BOOST_UBLAS_TYPE_CHECK
    BOOST_UBLAS_CHECK(ublas::detail::expression_type_check(ublas::prod(A, rhs), B), ublas::internal_logic());
#endif

    X.assign_temporary(rhs);
}

math.cpp に実一般係数行列の線形方程式系の解を計算する関数 solve を実装しました.math.cpp をプロジェクトに加え(他のソースコードと一緒にコンパイル),math.hpp を include することにより solve が使用可能になります.
solve は AX = B の解行列 X あるいは Ax = b の解ベクトル x を計算します.ただし,A は N 次正方行列,X と B は N x NRHS 行列(NRHS は 1 以上の任意の数),x と b は N 次元ベクトルを表します.


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

test.cpp

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


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

    ublas::matrix<double> A(3, 3);
    ublas::vector<double> b(3), x;

    A(0, 0) = 3; A(0, 1) = 2; A(0, 2) = -1;
    A(1, 0) = 2; A(1, 1) = -2; A(1, 2) = 4;
    A(2, 0) = -1; A(2, 1) = 0.5; A(2, 2) = -1;

    b(0) = 1; b(1) = -2; b(2) = 0;

    math::solve(A, b, x);

    // 以下結果の表示
    // 計算結果を http://en.wikipedia.org/wiki/System_of_linear_equations と比較してみましょう
    std::cout << x << std::endl;

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

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


次回は lapack を用いて同じ計算を実装する予定です.
それでは.