/*!
  \file    smdemo.cc
  \brief   demo for sparsematrix
  \date    4 May 2010
  \author  Clemens Pechstein

  Institute of Computational Mathematics
  Johannes Kepler University Linz, Austria

  Tutorials -- Numerical methods for elliptic PDEs
  summer semester 2010
*/
#include "vector.hh"
#include "sparsematrix.hh"

using namespace std;

int main ()
{
  cout << "vector demo" << endl
       << "===========" << endl << endl;
  
  Vector v;      // creates a vector of length 0
  v.resize (10); // resize
  
  v = 1; // set all entries to 0
  v[0] = 0;
  v[4] = -1.5;

  cout << "v = " << v << endl << endl;

  Vector w(10); // set size in constructor
  
  w = 0; // set all entries to 0
  w[1] = 4;
  w[6] = 7;
  w[9] = -3.3;
  
  cout << "w = " << w << endl << endl;

  cout << "||v|| = " << l2norm (v) << endl << endl;
  cout << "v . w = " << inner_product (v, w) << endl << endl;

  cout << endl;
  cout << "sparse shape demo" << endl
       << "=================" << endl << endl;

  SparseShape shape(10, 10); // pattern of a 10x10 matrix
  shape.addPosition (0, 0);  // make (0,0) non-zero element
  
  // create a trigagonal pattern
  for (int i=1; i<10; ++i)
    {
      shape.addPosition (i, i);
      shape.addPosition (i-1, i);
      shape.addPosition (i, i-1);
    }
  
  // add additional entries
  shape.addPosition (0, 6);
  shape.addPosition (2, 7);
  shape.addPosition (4, 8);
  shape.addPosition (5, 2);
  shape.addPosition (7, 3);
  shape.addPosition (8, 0);

  cout << "matrix pattern:" << endl;
  shape.print (cout);
  cout << endl;
  
  cout << endl;
  cout << "sparse matrix demo" << endl
       << "==================" << endl << endl;

  SparseMatrix K(shape); // create sparse matrix with the previously defined shape
  
  K = 0; // set all ("active") entries to zero;
  
  // set all diagonal entries to 2
  for (int i=0; i<10; ++i)
    K(i, i) = 2;
  
  // set all off-diagonal entries to -0.5
  for (int i=1; i<10; ++i)
    {
      K(i-1, i) = -0.5;
      K(i, i-1) = -0.5;
    }
  
  cout << "K has " << K.elems() << " non-zero/active entries" << endl << endl;

  cout << "K (short print): " << endl;
  K.printShort (cout);
  cout << endl;
  cout << "K (pretty print):" << endl;
  K.print (cout);
  cout << endl;
  // setting other entries
  K(0, 6) = -0.2;
  K(2, 7) = -0.1;
  K(4, 8) = -0.15;
  K(5, 2) = 0.01;
  K(7, 3) = -0.01;
  K(8, 0) = -0.13;
  
  cout << "K after additional settings:" << endl;
  K.print (cout);
  cout << endl;
    
  // What happens if accessing entries that are not in the pattern?

  // uncomment the following lines to get errors: this entry is not in the pattern
  //K(1, 9) = -1;
  //cout << "K(1, 9) = " << K(1, 9) << endl << endl;
  
  // for advanced users: the following would work
  //cout << "K(1, 9) = " << const_cast<const SparseMatrix&>(K)(1, 9) << endl << endl;
  
  return 0;
}