ExampleSparseMatrix2.cpp

/* This file is part of the BIAS library (Basic ImageAlgorithmS).

   Copyright (C) 2003-2009    (see file CONTACT for details)
   Multimediale Systeme der Informationsverarbeitung
   Institut fuer Informatik
   Christian-Albrechts-Universitaet Kiel


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

   BIAS is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU Lesser General Public License for more details.

   You should have received a copy of the GNU Lesser General Public License
   along with BIAS; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/



/**
  @example ExampleSparseMatrix2.cpp
  @relates SparseMatrix
  @brief   Example for sparse matrix with umfpack


  This examples demonstrates how the umfpack library can be used to solve
  sparse linear systems of equations.


  @see SparseMatrix.hh and umfpack userguide.pdf for more details on
       the sparse matrix methods

  @ingroup g_examples
  @author sedlazeck
  @date 06/2009
 */


#include <iostream>
#include <Base/Math/SparseMatrix.hh>
#include <MathAlgo/SVD.hh>
#include <Base/Math/Random.hh>

#ifdef BIAS_HAVE_UMFPACK
#  include <umfpack.h>
#endif

using namespace std;
using namespace BIAS;



int main()
{

  cout << "-- Testing umfpack - conversion and usage -- " << endl;

  cout << "does this get me anywhere  " << endl << flush;
  
  // matrix example from umfpack userguide
  BIAS::Matrix<double> A(5, 5, MatrixIdentity);
  A[0][0] =  2.0;
  A[0][1] =  3.0;
  A[1][0] =  3.0;
  A[1][1] =  0.0;
  A[1][2] =  4.0;
  A[1][4] =  6.0;
  A[2][1] = -1.0;
  A[2][2] = -3.0;
  A[2][3] =  2.0;
  A[3][2] =  1.0;
  A[3][3] =  0.0;
  A[4][1] =  4.0;
  A[4][2] =  2.0;
  A[4][4] =  1.0;
  BIAS::SparseMatrix S(A);

  cout << "sparse matrix S " << S << endl;

  Vector<double> b(5);
  b[0] =   8.0;
  b[1] =  45.0;
  b[2] =  -3.0;
  b[3] =   3.0;
  b[4] =  19.0;

  Vector<double> x(5);

  S.Solve(b, x);

  cout << "solved with sparse matrix (should be [1 2 3 4 5]) : " << x << endl << flush;

  SparseMatrix ST(S.GetColsNum(), S.GetRowsNum());
  S.Transpose(ST);
  cout << "ST " << ST << endl;

  // testing invert method on larger sparse matrix
  int m = 100;
  cout << "testing invert for larger matrix " << m << "x" << m << endl << flush;
  cout << "haaae? "<< endl; 
  Random random;
  Matrix<double> MLarge(m, m);
  double randElem;
  for(int i=0; i < m; i++){
    for(int j=0; j < m; j++){
      randElem = random.GetUniformDistributed(0.0, 1.0);
      if(randElem > 0.9){
        MLarge[i][j] = random.GetUniformDistributed(-100.0, 100.0);
        cout << "setting non-zero element " << endl;
      } else {
        MLarge[i][j] = 0.0;
      }
    }
  }
  // in case umfpack does not run, matrix needs to have last row and colum elem set
  MLarge[m-1][m-1] = random.GetUniformDistributed(-100.0, 100.0);
  
  cout << "new matrix is " << MLarge << endl << flush;
  
  cout << "building sparse matrix now... " << endl << flush;
  
  SparseMatrix SLarge(MLarge);
  SparseMatrix SLargeInv(SLarge.GetRowsNum(), SLarge.GetColsNum());

  cout << "maxCols " << SLarge.GetCols() << " actualCols " << SLarge.GetColsNum() 
       << " maxRows " << SLarge.GetRows() << " actual rows " << SLarge.GetRowsNum() << endl << flush;
 
  SLarge.Invert(SLargeInv);

  // reinit SLarge in case umfpack is not used
  SLarge = SparseMatrix(MLarge);
  
  Matrix<double> MLargeInv;
  SLarge.Invert(MLargeInv);


  cout << "Done with inverting, checking for id now... " << endl;

  SparseMatrix IdSparse(SLarge.GetRowsNum(), SLargeInv.GetColsNum());
  Matrix<double> IdSparseAsDense;
  SLarge.Multiply(SLargeInv, IdSparse);
  IdSparse.GetAsDense(IdSparseAsDense);

  cout << "Done getting sparse id as dense... " << endl;

  Matrix<double> IdDense;
  MLarge.Mult(MLargeInv, IdDense);


  double epsilon = 1e-8;
  bool res1 = IdSparseAsDense.IsIdentity(epsilon);
  bool res2 = IdDense.IsIdentity(epsilon);
  cout << "testing for identity - sparse : " << res1 << endl;
  cout << "testing for identity - dense  : " << res2 << endl;

  Matrix<double> ID(m, m, MatrixIdentity);
  cout << "checking norm sparse case " << Matrix<double>(IdSparseAsDense - ID).NormFrobenius() << endl;
  cout << "checking norm dense case " << Matrix<double>(IdDense - ID).NormFrobenius() << endl;


  SVD svd(MLarge);
  Matrix<double> MLargeInvSvd = svd.Invert(MLarge);
  Matrix<double> idMLarge;
  MLarge.Mult(MLargeInvSvd, idMLarge);
  cout << "checking norm in svd case svd res " << Matrix<double>(ID - idMLarge).NormFrobenius() << endl;

  if(!res1) exit(-1);
  if(!res2) exit(-1);

  exit(0);
}