ExampleSparseMatrix2.cpp

Example for sparse matrix with umfpackThis examples demonstrates how the umfpack library can be used to solve sparse linear systems of equations.

See Also

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

Author

sedlazeck

Date

06/2009

/* 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);
}