こんにちは uBLAS

C++ Boost uBLAS

大学のレポートでほげほげな数値計算してねー、とかいうのがあった。そんな課題がある度に毎回 C言語でうさんくさい行列計算ライブラリみたいの作ったりしてたけど、もうそんな無駄はやめよう。boost::numeric::ublas を使うことに。当該講義の教官は「うちの研究室では今 Boost のライブラリ使ってる」と言ってんで、読みなれてるでしょうし。C++ ってほとんど使ったこと無いけど、まあなんとかなんべ。
参考。

とりあえず単純なガウス消去（ピボット選択とかしないやつ）でも試してみる。

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

using namespace boost::numeric::ublas;

vector<double> gauss(matrix<double> a, vector<double> v)
{
  vector<double> result(3);
  for(int i=0;i<a.size1();++i){
    for(int j=0;j<a.size1();++j){
      if(i==j) continue;
      double mj = a(j,i)/a(i,i);
      for(int k=0;k<a.size2();++k){
        a(j,k) -= mj*a(i,k);
      }
      v(j) -= mj*v(i);
    }
  }
  for(int i=0;i<v.size();++i)
    result(i) = v(i)/a(i,i);
  return result;
}
int main()
{

  // 行列 (8 2 1)
  // 　　 (3 5 4)
  // 　　 (4 1 7)
  matrix<double> a( 3, 3 );
  a(0,0) = 8;
  a(0,1) = 2;
  a(0,2) = 1;
  a(1,0) = 3;
  a(1,1) = 5;
  a(1,2) = 4;
  a(2,0) = 4;
  a(2,1) = 1;
  a(2,2) = 7;

  // ベクトル (4 8 0)
  vector<double> v(3);
  v(0) = 4;
  v(1) = 8;
  v(2) = 0;

  std::cout << a << "x = " << v << std::endl;
  vector<double> x = gauss(a, v);
  std::cout << "x = " << x << std::endl;

  return 0;
}

出力

[3,3]((8,2,1),(3,5,4),(4,1,7))x = [3](4,8,0)
x = [3](0.0904977,1.79186,-0.307692)

えーと？

gosh> (define (f x1 x2 x3) (+ (* x1 .0904977) (* x2 1.79186) (* x3 -0.307692)))
f
gosh> (f 8 2 1)
4.000009599999999
gosh> (f 3 5 4)
8.000025100000002
gosh> (f 4 1 7)
6.799999999529405e-6

正しいみたいだ。よーし良かったー。

ではない。uBLAS を使ってレポートにとりくむのが目的だった。さて。LDV分解とかちまちま書くか。

test