Example for sparse matrix with umfpackThis examples demonstrates how the umfpack library can be used to solve sparse linear systems of equations.
#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;
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;
cout << "sparse matrix S " << S << endl;
b[0] = 8.0;
b[1] = 45.0;
b[2] = -3.0;
b[3] = 3.0;
b[4] = 19.0;
cout << "solved with sparse matrix (should be [1 2 3 4 5]) : " << x << endl << flush;
cout << "ST " << ST << endl;
int m = 100;
cout << "testing invert for larger matrix " << m << "x" << m << endl << flush;
cout << "haaae? "<< endl;
double randElem;
for(int i=0; i < m; i++){
for(int j=0; j < m; j++){
if(randElem > 0.9){
cout << "setting non-zero element " << endl;
} else {
MLarge[i][j] = 0.0;
}
}
}
cout << "new matrix is " << MLarge << endl << flush;
cout << "building sparse matrix now... " << endl << flush;
<<
" maxRows " << SLarge.
GetRows() <<
" actual rows " << SLarge.
GetRowsNum() << endl << flush;
cout << "Done with inverting, checking for id now... " << endl;
IdSparse.GetAsDense(IdSparseAsDense);
cout << "Done getting sparse id as dense... " << endl;
MLarge.
Mult(MLargeInv, IdDense);
double epsilon = 1e-8;
bool res1 = IdSparseAsDense.IsIdentity(epsilon);
cout << "testing for identity - sparse : " << res1 << endl;
cout << "testing for identity - dense : " << res2 << endl;
cout <<
"checking norm sparse case " <<
Matrix<double>(IdSparseAsDense - ID).NormFrobenius() << endl;
cout <<
"checking norm dense case " <<
Matrix<double>(IdDense - ID).NormFrobenius() << endl;
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);
}