dc/dee/ExampleSparseMatrix_8cpp-example.html

/*

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 ExampleSparseMatrix.cpp

    @relates SparseMatrix

    @brief Test implementation of sparse matrix operations.

    @see BIAS::SparseMatrix

    @ingroup g_examples

    @author esquivel

    @date 11/2008

*/


#include <iostream>

#include <Base/Math/Random.hh>

#include <Base/Math/SparseMatrix.hh>

#include <MathAlgo/SVD.hh>


using namespace std;

using namespace BIAS;


/** brief Helper function to get rid of small matrix entries. */

BIAS::Matrix<double> CropMatrix(const BIAS::Matrix<double> &M,

                                double epsilon = 1e-15)

{

  BIAS::Matrix<double> A(M);

  const unsigned int n = A.GetRows();

  const unsigned int m = A.GetCols();

  for (unsigned int i = 0; i < n; i++) {

    for (unsigned int j = 0; j < m; j++) {

      if (fabs(A[i][j]) < epsilon) A[i][j] = 0.0;

    }

  }

  return A;

}


int main()

{


  cout << "-- Testing base operations -- " << endl;


  BIAS::Matrix<double> A(6, 6, MatrixIdentity);

  A[2][1] = 17;

  A[3][5] = -0.1;

  A[4][4] = 0.0;

  BIAS::SparseMatrix S(A);


  cout << "- Dense is A = " << A << endl;

  cout << "- Sparse is S = " << endl;

  S.Print("S", cout);


  cout << "- Modifying S (consistency check)... " << flush;

  S.AddToDiagonal(2.0);

  cout << (S.CheckConsistency() ? "OK, " : "Not OK!, ") << flush;

  S.MultiplyDiagonalBy(0.5);

  cout << (S.CheckConsistency() ? "OK, " : "Not OK!, ") << flush;

  S.SetElement(1, 1, 0.0);

  cout << (S.CheckConsistency() ? "OK, " : "Not OK!, ") << flush;

  S.SetElement(1, 4, 1.0);

  cout << (S.CheckConsistency() ? "OK" : "Not OK!")

       << endl << flush;

//

//  S.SetElement(6,2,1.0);

//  S.SetElement(6,2,2.0);

//  cout << (S.CheckConsistency() ? "OK" : "Not OK!")

//  << endl << flush;


  cout << "sparse S after setting single elements is " << S << endl << flush;


  SparseMatrix SAddTest(2,2);

  SAddTest.SetElement(0,0,1.0);

  SAddTest.SetElement(1,1,1.0);

  cout << "s before adding " << SAddTest << endl;

  Matrix<double> addBottom(2,3,MatrixIdentity);

  addBottom[0][2] = 2;

  Matrix<double> addRight(3,2,MatrixIdentity);

  addRight[2][0] = 3;

  SparseMatrix addBottomSparse(addBottom);

  SparseMatrix addRightSparse(addRight);


  cout << "adding to bottom " << addBottomSparse  << addBottomSparse.GetRows() << " " << addBottomSparse.GetCols() << endl;

  cout << "adding to right " << addRightSparse << addRightSparse.GetRows() << " " << addRightSparse.GetCols() <<endl;

  SAddTest.AppendMatrixBottom(addBottomSparse);

  cout << "s after adding to bottom " << SAddTest << " size " << SAddTest.GetRows() << " " << SAddTest.GetCols() << endl;

  SAddTest.AppendMatrixRight(addRightSparse);

  cout << "s after adding to right " << SAddTest << " size " << SAddTest.GetRows() << " " << SAddTest.GetCols() << endl;


  cout << (SAddTest.CheckConsistency() ? "OK" : "Not OK!")

  << endl << flush;


  BIAS::Matrix<double> B;

  S.GetAsDense(B);


  cout << "- Modified sparse is S = " << endl;

  S.Print("S", cout);

  cout << "- Conversion to dense is B = " << B << endl;


  BIAS::Matrix<double> BT, BTB;

  BT = B.Transpose();

  BTB = BT * B;

  BIAS::SparseMatrix ST(S.GetColsNum(), S.GetRowsNum()), STS(ST.GetRowsNum(), S.GetColsNum());

  S.Transpose(ST);

  ST.Multiply(S, STS);


  cout << "- Sparse multiplication S'S = " << endl;

  STS.Print("S'S", cout);

  cout << "- Dense multiplication B'B = " << BTB << endl;


  double m;

  unsigned int i, j;

  m = STS.GetMaxDiagonalElement(i);

  cout << "- Maximal diagonal element of S'S is [" << i << "][" << i

       << "] = " << m << endl;

  j = 4;

  m = S.GetMaxRowElement(j, i);

  cout << "- Maximal element of row " << j << " in S is [" << j << "][" << i

       << "] = " << m << endl;

  m = ST.GetMaxColumnElement(j, i);

  cout << "- Maximal element of column " << j << " of S' is [" << i

       << "][" << j << "] = " << m << endl;


  cout << "- Checking consistency of S... ";

  cout << (S.CheckConsistency() ? "OK!" : "failed!") << endl;

  cout << "- Checking consistency of S'... ";

  cout << (ST.CheckConsistency() ? "OK!" : "failed!") << endl;

  cout << "- Checking consistency of S'S... ";

  cout << (STS.CheckConsistency() ? "OK!" : "failed!") << endl;


  cout << endl << "-- Testing matrix inversion -- " << endl << endl;


  BIAS::Matrix<double> M(6, 6, MatrixZero);

  M[0][0] =  2.0;

  M[1][1] =  4.0;

  M[2][2] = -1.0;

  M[3][3] = -3.0;

  M[4][4] =  0.5;

  M[5][5] =  1.0;

  M[3][4] =  1.5;

  M[2][4] =  0.1;

  M[3][1] =  2.5;

  BIAS::SparseMatrix N(M);

  BIAS::SparseMatrix Nold(M);

  N.SetElement(5, 5, 1.0);

  cout << "- Dense is M = " << M << endl;

  cout << "- Sparse is N = " << endl;

  N.Print("N", cout);

  cout << "- Checking consistency of N... ";

  cout << (N.CheckConsistency() ? "OK!" : "failed!") << endl;


  BIAS::SVD svd;

  BIAS::Matrix<double> Minv = svd.Invert(M);

  cout << "SVD result of inverse is M^-1 = " << CropMatrix(Minv) << endl;


  BIAS::SparseMatrix Ninv(N.GetRowsNum(), N.GetColsNum());

  N.Invert(Ninv);


  cout << "- Sparse after inversion by Gauss-Jordan is N = " << endl;

  N.Print("N", cout);

  cout << "- Inverse is N^-1 = " << endl;

  Ninv.Print("N^-1", cout);

  cout << "- Checking consistency of N... ";

  cout << (N.CheckConsistency() ? "OK!" : "failed!") << endl;

  cout << "- Checking consistency of N^-1... ";

  cout << (Ninv.CheckConsistency() ? "OK!" : "failed!") << endl;


  cout << "- M * M^-1 = " << CropMatrix(M*Minv) << endl;


  BIAS::SparseMatrix NNinv(Nold.GetRowsNum(), Ninv.GetColsNum());

  Nold.Multiply(Ninv, NNinv);

  cout << "- N * N^-1 = " << endl;

  NNinv.Print("N*N^-1", cout);

  cout << "- Checking consistency of N*Ninv... ";

  cout << (NNinv.CheckConsistency() ? "OK!" : "failed!") << endl;


  cout << endl << "-- Testing linear equation systems -- " << endl << endl;


  BIAS::SparseMatrix L(M);

  BIAS::SparseMatrix Lold(M);

  BIAS::Vector<double> b(6), x(6), b2(6), x2(6);

  b[0] =  1.0;

  b[1] =  2.0;

  b[2] =  3.0;

  b[3] = -4.0;

  b[4] =  5.0;

  b[5] = -6.0;


  cout << "- Solving linear equation system L*x = b via dense SVD and sparse "

       << "matrix Gauss-Jordan approach..." << endl;

  svd.Solve(M, b, x);

  L.Solve(b, x2);

  cout << "- Sparse LES matrix after Gauss-Jordan is L = " << endl;

  L.Print("L", cout);

  cout << "- Checking consistency of L... ";

  cout << (L.CheckConsistency() ? "OK!" : "failed!") << endl;


  cout << "- Result of SVD is x = " << x << endl;

  cout << "- Result of sparse Gauss-Jordan is x = " << x2 << endl;

  Lold.Multiply(x2, b2);

  cout << "- Evaluation of |L*x - b| yields "

       << BIAS::Vector<double>(b-b2).NormL2()

       << " (compare with dense SVD: "

       << BIAS::Vector<double>(b-(M*x)).NormL2()

       << ")" << endl << endl;


  return 0;

}