ublasのbanded_matrixとlu_factorizeで、帯行列の逆行列とベクトルの積を線形時間で

帯行列の逆行列は一般には密行列なので、帯行列の逆行列の計算は。しかし、帯行列のLU分解はで出来るので、帯行列とベクトルに対してはで計算できる。これを実装。

/*! g++ main.cpp -g -DNDEBUG
 */

#include <boost/numeric/ublas/banded.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/operation.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>

#include <boost/timer.hpp>

#include <vector>
#include <iterator>

using namespace boost::numeric;
using namespace std;


void makeBanded(ublas::banded_matrix<double>& a,
		const int n, const int w){
  a.resize(n, n, w, w);
  a.clear();
  for(int i = 0; i < n; ++i){
    for(int j = max(0, i - w); j < i; ++j){
      a(i, j) = 1;
    }
    for(int j = i + 1; j < min(n, i + w + 1); ++j){
      a(i, j) = 2;
    }
    a(i, i) = w * 5;
  }
}

template<class V>
void makeVec(ublas::vector_expression<V>& v, size_t i){
  //copy(v().begin(), v().end(), 0); //compile error : the cause is unclear.
  for(size_t j = 0; j < v().size(); j++){
    v()(j) = 0;
  }
  v()(i) = 1;
}

template<class V>
void invertProd(ublas::banded_matrix<double>& a,
		ublas::vector_expression<V>& v){
  ublas::permutation_matrix<> pm(a.size1());
  ublas::lu_factorize(a, pm);
  ublas::lu_substitute(a, pm, v());
}




int main(int argc, char* argv[]){
  if(argc != 3){
    cout << "usage : a.out dim iter" << endl;
    return 1;
  }

  const size_t dim = atoi(argv[1]);
  const size_t iter = atoi(argv[2]);
  const size_t wid = 1;

  boost::timer timer;

  ublas::banded_matrix<double> a;
  makeBanded(a, dim, wid);

  ublas::vector<double> v(a.size1());
  makeVec(v, 0);

  ublas::banded_matrix<double> b(a);
  ublas::vector<double> ret(v);

  timer.restart();
  invertProd(b, ret);
  cout << timer.elapsed() << endl;

#ifndef NDEBUG
  ublas::matrix<double> aInv(a.size1(), a.size2());
  for(size_t i = 0; i < aInv.size1(); ++i){
    ublas::banded_matrix<double> b(a);
    ublas::matrix_column<ublas::matrix<double> > aInvCol(aInv, i);
    makeVec(aInvCol, i);
    invertProd(b, aInvCol);
  }

  cout << a << endl
       << prod(a, aInv) << endl;
#endif

  return 0;
}

#ifndef NDEBUGの中は、invertProdが正しく逆行列との積を計算していることの確認。の結果をaInvに入れるている。これで、の第列を取り出して、一つの行列にまとめた上でとの積を計算して、単位行列になることを確認。
関数の引数をvector_expression型にすると、vectorそのものでも、matrix_columnのようなものでも受け取れるということは分かったんだけど、イテレータの部分がコンパイルエラーになったりして、イマイチ使い方が良く判らん。こういうのは、どこで勉強すれば良いのか。