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

今回は前回の連立方程式の計算ã‚’ lapack を用いて実装します．
lapack には dgesv という関数があり，これが内部で係数行列をLU分解して解を計算してくれます．私たちがするべきことといえば，いつものように bindings.lapack.gesv に渡す column_major な行列を準備することくらいでしょう．以下に私の実装例を示します．

math.hpp

#ifndef MATH_HPP_20081106
#define MATH_HPP_20081106

#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/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/vector.hpp>
#define BOOST_NUMERIC_BINDINGS_USE_CLAPACK
#include <boost/numeric/bindings/lapack/gesv.hpp>
#include <boost/numeric/bindings/traits/ublas_matrix.hpp>
#undef  BOOST_NUMERIC_BINDINGS_USE_CLAPACK


void math::solve(const boost::numeric::ublas::matrix<double>& A, const boost::numeric::ublas::vector<double>& b, boost::numeric::ublas::vector<double>& x)
{
    BOOST_UBLAS_CHECK(A.size1() == A.size2(), boost::numeric::ublas::external_logic());
    BOOST_UBLAS_CHECK(A.size1() == b.size(), boost::numeric::ublas::external_logic());

    boost::numeric::ublas::matrix<double> B(b.size(), 1), X;

    boost::numeric::ublas::column(B, 0).assign(b);
    solve(A, B, X);

    x = boost::numeric::ublas::column(X, 0);
}


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, ublas::column_major> CA(A), CX(B);
    int                                        info;

    info = boost::numeric::bindings::lapack::gesv(CA, CX);
    BOOST_UBLAS_CHECK(info == 0, ublas::internal_logic());

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

    X = CX;
}

math.cpp に実一般係数行列の線形方程式系の解を計算する関数 solve を実装しました．デバッグ時には実際に計算した X を用いて，AX = B になることを確認するコードを追加しています．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 次元ベクトルを表します．
NRHS が 1 のとき，X と B は N x 1 の行列になりますが，これは N 次元ベクトルに他なりません（NRHS が 1 のときは，AX = B という式はベクトルを用いて Ax = b になるのです）．uBLAS ではベクトルを表す型もちゃんと用意されているため，引数にベクトルを取る関数(Ax = b を解く関数)も実装しました．コードを見ると分かりますが，この関数は AX = B を解く関数を利用しているだけです．

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

test.cpp

#include "math.hpp"
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.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> 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)

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

次回は対称係数行列の線形方程式系の解を計算する関数を実装する予定です．
それでは．