banded_matrixの積はちゃんと線形時間だった件

http://d.hatena.ne.jp/hotoku/20101201/1291219911の疑問が解決した。簡単なことだ。NDEBUGマクロを定義すれば良かったのだ。ublasは、NDEBUGが無い場合はデバッグビルドと看做してサイズチェック等のテストコードを埋め込む。なので、最大の速度を求めるならNDEBUGマクロをプリプロセッサ等で定義しないとダメ。ちなみにプリプロセッサはg++なら-DNDEBUG。
一応、lu_factorize(LU分解)とlu_substitute(後退代入と前進消去)関数でも線形時間であることを確かめた。
テスト用の実行ファイルソース。

/*! 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/timer.hpp>

#include <vector>
#include <iterator>

using namespace boost::numeric;
using namespace std;


class Check{
public:
  Check(size_t d, size_t i) : dim(d), iter(i), wid(1){
    check();
  }

private:
  const size_t dim, iter, wid;
  boost::timer timer;
  vector<double> time;

  void check(){
    band_vec();
    band_band();
    band_band_dense();
    lu();
    print();
  }

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

    ublas::vector<double> v1(a.size2());
    makeVec(v1);

    ublas::vector<double> v2(a.size1());

    start();
    for(size_t i = 0; i < iter; ++i){
      noalias(v2) = prod(a, v1);
    }
    end();
  }

  void band_band(){
    ublas::banded_matrix<double> a, b;
    makeBanded(a, dim, wid);
    makeBanded(b, dim, wid);

    ublas::banded_matrix<double> c(a.size1(),
				   b.size2(),
				   a.lower() + b.lower(),
				   a.upper() + b.upper());

    start();
    for(size_t i = 0; i < iter; ++i){
      noalias(c) = prod(a, b);
    }
    end();
  }

  void band_band_dense(){
    ublas::banded_matrix<double> a, b;
    makeBanded(a, dim, wid);
    makeBanded(b, dim, wid);

    ublas::matrix<double> c(a.size1(), b.size2());

    start();
    for(size_t i = 0; i < iter; ++i){
      noalias(c) = prod(a, b);
    }
    end();    
  }

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

    ublas::banded_matrix<double> b(a.size1(), a.size2(),
				   a.lower(), a.upper());
    b.assign(a);

    ublas::permutation_matrix<> pm(a.size1());

    start();
    for(size_t i = 0; i < iter; ++i){
      ublas::lu_factorize(b, pm);
    }
    end();

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

    start();
    for(size_t i = 0; i < iter; ++i){
      lu_substitute(b, pm, v);
    }
    end();
  }


  void start(){
    timer.restart();
  }

  void end(){
    time.push_back(timer.elapsed());
  }

  void print(){
    copy(time.begin(), time.end(), ostream_iterator<double>(cout, "\t"));
    cout << endl;
  }


  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 < min(n, i + w + 1); ++j){
	a(i, j) = 1;
      }
      a(i, i) += 3*w;
    }
  }

  void makeVec(ublas::vector<double>& v){
    for(size_t i = 0; i < v.size(); ++i){
      v(i) = i;
    }
  }


};


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

  Check(atoi(argv[1]),
	atoi(argv[2]));

  return 0;
}

帯行列と帯行列の積を受け取る変数を普通の行列にしちゃうとかかるので注意。
テストを好きなだけ繰り返す

#!/usr/bin/env python

from subprocess import PIPE, Popen

dat = open("dat.txt", "w")

step = 200
n = step
m = 1
nMax = 20000
while n <= nMax:
    proc = Popen(["./a.out", str(n), str(m)], stdout=PIPE, stderr=PIPE)
    (o, e) = proc.communicate()
    print n
    dat.write("%i\t%s\t%s" % (n, o, e))

    n += step

dat.close()    

テスト結果を解析

dat <- read.table("dat.txt")

for(i in 2:ncol(dat)){
  print(lm(log(dat[,i]) ~ log(dat[,1])))
}