Test implementation of sparse matrix operations.
#include <iostream>
#include <Base/Math/Random.hh>
#include <Base/Math/SparseMatrix.hh>
#include <MathAlgo/SVD.hh>
using namespace std;
using namespace BIAS;
double epsilon = 1e-15)
{
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;
A[2][1] = 17;
A[3][5] = -0.1;
A[4][4] = 0.0;
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;
cout << "sparse S after setting single elements is " << S << endl << flush;
SAddTest.SetElement(0,0,1.0);
SAddTest.SetElement(1,1,1.0);
cout << "s before adding " << SAddTest << endl;
addBottom[0][2] = 2;
addRight[2][0] = 3;
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;
S.GetAsDense(B);
cout << "- Modified sparse is S = " << endl;
cout << "- Conversion to dense is B = " << B << endl;
BTB = BT * B;
BIAS::SparseMatrix ST(S.GetColsNum(), S.GetRowsNum()), STS(ST.GetRowsNum(), S.GetColsNum());
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;
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;
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;
cout << "SVD result of inverse is M^-1 = " << CropMatrix(Minv) << endl;
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;
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;
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 "
<< " (compare with dense SVD: "
<< ")" << endl << endl;
return 0;
}