Re: [ublas] possible bug (probably in diagonal_matrix)

From: Gunter Winkler (guwi17_at_[hidden])
Date: 2005-06-24 07:01:38

On Friday 24 June 2005 10:47, Gunter Winkler wrote:
> > >         diagonal_matrix<double> z(sz), x(sz);
> > >         matrix<double> y(sz, sz);
> > >         matrix<double> B = prod(z, x - y);

Just another remark:

I would strongly discourage the use of the expression
matrix<double> B = prod(z, x - y);

because the computation is _very_ slow. (The banded iterators produce a big 
overhead on a dense matrix.)

gcc-3.3 times:
$ ./diag 1000
fill: 0.03
temp 0.39
direct 16.62
sum 0.43
diff = 0 0

gcc-3.4 times:
$ ./diag 1000
fill: 0.02
temp 0.25
direct 10.67
sum 0.33
diff = 0 0

the program:
#include <boost/timer.hpp>

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/banded.hpp>
#include <boost/numeric/ublas/io.hpp>

int main(int argc, char *argv[]) 
{
  using namespace boost::numeric::ublas;

  size_t sz = 1000;
  if (argc > 1)
    sz = ::atoi(argv [1]);

  {
    boost::timer t;

    t.restart();
    diagonal_matrix<double> z(sz), x(sz);
    matrix<double> y(sz, sz);
    for (size_t i = 0; i < sz; ++i) {
      z(i,i) = i;
      x(i,i) = 2*i;
      for (size_t j = 0; j < sz; ++j) {
        y(i,j) = i + j;
      }
    }
    std::cout << "fill: " << t.elapsed() << "\n";

    t.restart();
    matrix<double> x_y = x - y;
    matrix<double> A = prod(z, x_y);
    std::cout << "temp " << t.elapsed() << "\n";

    t.restart();
    matrix<double> B = prod(z, x - y);
    std::cout << "direct " << t.elapsed() << "\n";

    t.restart();
    matrix<double> C = prod(z, x) - prod(z, y);
    std::cout << "sum " << t.elapsed() << "\n";

    std::cout << "diff = " << norm_inf(A-B) << " " << norm_inf(A-C) << "\n";
  }

  return 0;
}

mfg
Gunter