ublastest.cpp

Go to the documentation of this file.

/* main.cpp  -  main routines */
/***************************************************************************
    begin                : Don Mär 14 10:19:20 MET 2002
    copyright            : (C) 2002 by Gunter Winkler
    email                : guwi17@gmx.de
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

/*
 * $Log$
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <iostream>
#include <iomanip>
#include <sstream>
#include <fstream>

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

#include "itl/interface/boost.h"
#include "itl/preconditioner/ilu_new.h"

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/operation.hpp>


/* private includes */

#ifndef NDEBUG
volatile const static char _main_RCSID[] = "$Id$";
#endif

using std::cout;
using std::cin;
using std::endl;
using std::ofstream;

using std::min;
using std::max;

using namespace boost::numeric::ublas;

typedef sparse_matrix<double> uBLASSparseMatrix;
typedef vector<double>            uBLASVector;

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

  // type       of stiffnes matrix
  typedef sparse_matrix<double>   MY_STIMA;
  typedef matrix_row< MY_STIMA >  MY_ROW;
  typedef vector<double>          MY_VEC;
  timeval t[7];

  int N = 100;
  //  cin >> N;

  MY_STIMA M(4,4,16);

  for (int i=0; i<M.size1(); ++i) M(i,i)=1.0;
  M(0,0)=3;
  M(0,1)=-1;
  M(0,2)=-2;
  M(1,0)=-1;
  M(1,1)=4;
  M(1,2)=-2;
  M(2,0)=-2;
  M(2,1)=-2;
  M(2,2)=9;

  itl::ILU<MY_STIMA> precond(M);

  MY_VEC xx(M.size2());
  MY_VEC yy(M.size2());

  yy(0) = 0; yy(1) = 1; yy(2) = 5;

  precond().solve(yy,xx);

  cout << "x= " << xx << "\n";

//   for (MY_STIMA::iterator1 i=M.begin1(); i!=M.end1(); ++i) {
//      cout << i.index1() << " " << i.index2() << " " << (*i) << "\n";
//   }

  MY_ROW r(M, 2);

  for (MY_ROW::iterator i=r.begin(); i!=r.end(); ++i) {
        cout << i.index() << " " << (*i) << "\n";
  }

  MY_STIMA::iterator1 r_it1 = M.begin1();
  MY_STIMA::iterator2 c_it1 = M.begin2();
  while ( (r_it1 != M.end1()) && (c_it1 != M.end2()) ) {
        MY_STIMA::iterator1  r_it2 = c_it1.begin();
        MY_STIMA::iterator2  c_it2 = r_it1.begin();
        MY_STIMA::iterator2  c;

        c = c_it2;
        while ( (r_it2 != c_it1.end()) && (c != r_it1.end()) ) {
          while (c.index2() < r_it2.index1()) { ++c; }
          while (c.index2() > r_it2.index1()) { ++r_it2; }
          if (c.index2() == r_it2.index1()) {
                cout << r_it2.index1() << ", " << c.index2() << ": " << (*c) << " += " << (*r_it2);
                (*c) += (*r_it2);
                cout << " --> " << (*c) << " ? \n";
                ++ r_it2;
                ++ c;
          }
        }
        ++ r_it1;
        ++ c_it1;
  }

  cout << "M = " << M << endl;

  cout << "Benchmark for linear algebra in uBLAS" << endl;
  cout << "--------------------------------------" << endl;

  gettimeofday( &t[0], NULL );

  cout << "Creating a vector with " << N << " elements" << endl;
  uBLASVector x(N);

  gettimeofday( &t[1], NULL );

  cout << "Creating another vector with " << N << " elements" << endl;
  uBLASVector y(N);

  gettimeofday( &t[2], NULL );

  cout << "Creating a " << N << " x " << N << " matrix" << endl;
  uBLASSparseMatrix A(N, N);

  gettimeofday( &t[3], NULL );

  cout << "Assembling" << endl;
  for (int i = 0; i < N; i++) {
    A(i,i) += 2.0;

    if ( i > 0 )
      A(i, i - 1) += -1.0;

    if ( i < (N-1) )
      A(i, i + 1) += 1.0;
  }

  gettimeofday( &t[4], NULL );

  cout << "Multiplying: y = A*x" << endl;
  axpy_prod(A, x, y);

  gettimeofday( &t[5], NULL );

  double res = 0;
  for (int i = 0; i<N; ++i) res += y[i];
  cout << res << endl;

  for (int i = 0; i<5; ++i) {
        double t1=(double(t[i].tv_sec)+t[i].tv_usec/1000000.0);
        double t2=(double(t[i+1].tv_sec)+t[i+1].tv_usec/1000000.0);
        double time=(t2-t1);
        cout << "time " << i << ": " << time << "\n";
  }
  cout << endl;

  return EXIT_SUCCESS;
};

---