mtl_trisolve.cpp

Go to the documentation of this file.

#include <sys/time.h>
#include <unistd.h>

#include <iostream>

#include "mtl/mtl.h"
#include "mtl/matrix.h"
#include "mtl/mtl_algo.h"

using std::cout;
using std::endl;

using namespace mtl;

/*
 * LU for matrix A
 *
 * for i = 2:n
 *   for k = 1:i-1
 *     A(i,k) = A(i,k) / A(k,k);
 *     for j = k+1:n
 *       A(i,j) = A(i,j) - A(i,k)*A(k,j)
 *     next j
 *   next k
 * next i
 *
 * we assume: all diagonal elements are nonzero
 */
    
template <class Matrix>
void do_lu(Matrix& LU)
{
  typedef  typename Matrix::size_type   size_type;

  typename Matrix::iterator r = LU.begin();
  typename Matrix::iterator r_end = LU.end();

  ++r;
  while (r != r_end) {
        size_type  i = r.index();
        typename Matrix::iterator rk = LU.begin();
        typename Matrix::iterator rk_end = LU.end();
        typename Matrix::OneD::iterator c_start = (*r).begin();
        typename Matrix::OneD::iterator c_end = (*r).end();
        typename Matrix::OneD::iterator c;
        typename Matrix::OneD::iterator kk;
        while (c_start.index() < i) {
          size_type  k = c_start.index();
          while (rk.index() < k) ++rk;
          c = c_start;
          kk = (*rk).begin();
          while (kk.index() < k) ++kk; // better use find2(...) ?
          (*c) /= (*kk);
          ++c; ++kk;
          while (kk != (*rk).end()) {
                // insert elements if necessary
                LU(r.index(),kk.index()) -= (*c_start) * (*kk);
                ++kk;
          }
          ++c_start;
        }
        ++r;
  }
}

template <class Matrix>
void spy(std::ostream& os, const Matrix& A)
{
  const double eps = 1e-10;
  typename Matrix::size_type i,j,m,n;
  m = A.nrows();
  n = A.ncols();
  for (i=0; i<m; ++i) {
        for (j=0; j<n; ++j) {
          if (A(i,j) < -eps) {
                os << "-";
          } else if (A(i,j) > eps) {
                os << "+";
          } else /* if A(i,j) == 0 */ {
                os << " ";
          }
        }
        os << "\n";
  }
  os << "\n";
}

int main(int argc, char *argv[])
{

  typedef mtl::matrix< double, 
        rectangle<>, 
        mtl::array< mtl::compressed<> >, 
        row_major >::type  MY_STIMA;

  typedef mtl::dense1D<double>         MY_VEC;

  const int COUNT= 2;

  const int size = 20;
  const int sub1 = 1;
  const int sub2 = 10;
  const int sup1 = 1;
  const int sup2 = 10;

  timeval t[7];

  // MY_STIMA M(size,size,(sup2+1+sub2)*size);
  MY_STIMA M(size,size);
  MY_VEC   x(size);
  MY_VEC   y(size);
  set(y,0.0);

  for (int i=sup2; i<size; ++i)  M(i-sup2,i     ) =-1.0;
  for (int i=sup1; i<size; ++i)  M(i-sup1,i     ) =-1.0;
  for (int i=0;    i<size; ++i)  M(i     ,i     ) = 4.0;
  for (int i=sub1; i<size; ++i)  M(i     ,i-sub1) =-1.0;
  for (int i=sub2; i<size; ++i)  M(i     ,i-sub2) =-1.0;

  //cout << "input M = " << M << "\n";
  if (size < 72) {
        cout << "matrix spy of input matrix\n";
        spy(cout, M);
  }

  for (int i=0; i<x.size(); ++i)  x[i] = 1.0;

  gettimeofday( &t[0], NULL );
  for (int i=0; i<COUNT; ++i) {
        mult(M, x, y);
  }
  gettimeofday( &t[1], NULL );

  /* uncomment LU if you want to check the results, but my LU is slow ... */
  do_lu(M);
  gettimeofday( &t[2], NULL );

  if (size < 72) {
        cout << "matrix spy after LU decomposition\n";
        spy(cout, M);
  }
  // cout << "output M = "; print_all_matrix(M);
  
  cout << "start solver\n";

  // note: mtl-trisolve fails if there are elements in the wrong half ...
  gettimeofday( &t[3], NULL );
  tri_solve(triangle_view<MY_STIMA, unit_lower>::type (M), y);
  print_vector(y);
  gettimeofday( &t[4], NULL );
  tri_solve(triangle_view<MY_STIMA, upper>::type (M), y);
  print_vector(y);
  gettimeofday( &t[5], NULL );
  
  add(scaled(y, -1.0), x);
  cout << "L2 error : " << (two_norm(x)) << "\n";

  cout << "### running times " << size << endl;
  double t1,t2;
  t1=(double(t[0].tv_sec)+t[0].tv_usec/1000000.0);
  t2=(double(t[1].tv_sec)+t[1].tv_usec/1000000.0);
  cout << "### axpy_prod   : " << (t2-t1)/COUNT << " s " << endl;
  t1=(double(t[1].tv_sec)+t[1].tv_usec/1000000.0);
  t2=(double(t[2].tv_sec)+t[2].tv_usec/1000000.0);
  cout << "### compute LU  : " << (t2-t1) << " s " << endl;
  t1=(double(t[3].tv_sec)+t[3].tv_usec/1000000.0);
  t2=(double(t[4].tv_sec)+t[4].tv_usec/1000000.0);
  cout << "### solving unit lower system : " << (t2-t1) << " s " << endl;
  t1=(double(t[4].tv_sec)+t[4].tv_usec/1000000.0);
  t2=(double(t[5].tv_sec)+t[5].tv_usec/1000000.0);
  cout << "### solving upper system      : " << (t2-t1) << " s " << endl;

  return EXIT_SUCCESS;
};

---