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 あたり参考に適当に書いたんだけど。詳しい意味は知らない。たぶんきっと「どんな型でもかかってこいや!」って意味だと思う。