uBLAS べんりだ。

C++ uBLAS

std::cin から行列、ベクトルへの入力も用意されてる。

% cat helloublas.cpp
#include <iostream>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>

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

  matrix<int> a;
  vector<int> v;
  std::cin >> a >> v;
  std::cout << "a = " << a << std::endl << "v = " << v << std::endl;

  return 0;
}
% g++ helloublas.cpp
% ./a.out
[2,3]((1,2,3),(4,5,6))
[5](1,2,3,4,5)
a = [2,3]((1,2,3),(4,5,6))
v = [5](1,2,3,4,5)

C++ の stream とかどうなんだろうなあと思ってたけど、この辺は便利。この辺までやらんと stdio のが便利な気するけど。

ついでにこの辺使ったガウス消去とかも書いた。

% cat in_gauss_elim1
[3,3]((8,2,1),(3,5,4),(4,1,7))
[3](4,8,0)
% ./gauss_elim <in_gauss_elim1
[3,3]((8,2,1),(3,5,4),(4,1,7))x = [3](4,8,0)
x = [3](0.0904977,1.79186,-0.307692)

ちゃんとできてる。

以下コード。

#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;

template<class M, class V> void gauss(M& a, V& v)
{
  for(size_t i=0;i<a.size1();++i){
    for(size_t j=0;j<a.size1();++j){
      if(i==j) continue;
      double mj = a(j,i)/a(i,i);
      for(size_t k=0;k<a.size2();++k){
        a(j,k) -= mj*a(i,k);
      }
      v(j) -= mj*v(i);
    }
  }
  for(size_t i=0;i<v.size();++i)
    v(i) /= a(i,i);
}
int main()
{
  matrix<double> a;
  vector<double> v;
  std::cin >> a >> v;

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

  return 0;
}

-Wall 付きでコンパイルしたら

for(int i=0;i<a.size();++i){
  // ...
}

の部分で「符号付きと符号無しの数比較してんじゃねえよタコ！」とか言われたので size_t とか使ってみた。でもこれだと降順で処理したいときについ

for(size_t i=a.size();i>=0;--i){
  // ...
}

みたいな無限ループを書いてしまいそうで怖い。というか書いた。スマートな解決法てなんだ。

for(int i=0;i<(int)a.size();++i){
  // ...
}

ですか。たぶん違う気がする。

あと template ってなんでしょうね。matrix<double>& aと書けなかったから matrix.hpp あたり参考に適当に書いたんだけど。詳しい意味は知らない。たぶんきっと「どんな型でもかかってこいや！」って意味だと思う。

test