ublasのlu_factorizeはbanded_matrixには使えないのか?
http://d.hatena.ne.jp/hotoku/20110121/1295538662なんて書いて喜んだのも束の間、ublasで帯行列とLU分解の組み合わせは一筋縄では行かぬかもしれない。
LU分解関数 lu_factorizeの中のpivotingの部分で、matrix_rowクラスのswap関数を呼ぶんだけど、こいつが帯行列相手にはきちんと機能しないらしい。
以下テスト。
/*! g++ main.cpp -g */ #include <boost/numeric/ublas/banded.hpp> #include <boost/numeric/ublas/matrix_proxy.hpp> #include <boost/numeric/ublas/io.hpp> #include <iostream> using namespace std; namespace ublas = boost::numeric::ublas; int main(int argc, char* argv[]){ typedef boost::numeric::ublas::banded_matrix<double> ublas_banded; const int n = 4; const int w = 2; ublas_banded mat(n, n, w, w); for(int i = 0; i < n; ++i){ for(int j = max(0, i - w); j < min(n, i + w + 1); ++j){ mat(i, j) = i; } } cout << "before : " << mat << endl; ublas::matrix_row<ublas_banded> mr(ublas::row(mat, n - 2)); ublas::row(mat, n - 1).swap(mr); cout << "after : " << mat << endl; return 0; }
出力はこれ。
bash-3.2$ ./a.out before : [4,4]((0,0,0,0),(1,1,1,1),(2,2,2,2),(0,3,3,3)) after : [4,4]((0,0,0,0),(1,1,1,1),(2,3,3,3),(0,2,2,2))
2行目と3行目(0始まりindexで)を交換したつもりなんだけど、交換出来てない。ただ、今回のデバッグ作業で気が付いたんだけど、pivoting無しバージョンのlu_factorize関数もあるらしい。取り敢えず、正定値行列扱う分にはそれ使っとこう。
日記見ると、今回のデバッグだけで一週間くらい掛かっているじゃないか。チキショウ、時間の無駄だ。自分には時間が無いというのに。