#include <bias_config.h>
#include <MathAlgo/SVD.hh>
#include <MathAlgo/SVD3x3.hh>
#include <Base/Math/Random.hh>
#include <Base/Math/Matrix3x3.hh>
#include <Base/Debug/TimeMeasure.hh>
using namespace BIAS;
using namespace std;
#define MAX_PRINT_SIZE 8
int main(int argc, char **argv)
{
int size = 3;
if (argc > 1) {
size = atoi(argv[1]);
}
bool useOpenCV = false;
if (argc > 2)
useOpenCV = (atoi(argv[2]) != 0);
if (size < 1) {
cout << "Usage : ExampleSVD [<matrix size>] [<use OpenCV?>]" << endl << endl;
return 0;
}
if (useOpenCV) {
#ifdef BIAS_HAVE_OPENCV
cout << endl << "Using OpenCV singular value decomposition" << endl;
#else
cout << "BIAS does not support OpenCV, using LAPACK SVD instead" << endl;
useOpenCV = false;
#endif
} else {
cout << endl << "Using LAPACK singular value decomposition" << endl;
}
double errorThreshold = 1e-6 * size * size;
cout << endl << "1.) Compute SVD for random " << size << "x" << size << " matrix" << endl;
for (int i = 0; i < size; i++)
for (int j = 0; j < size; j++)
if (size <= MAX_PRINT_SIZE)
int numRepeats = 100;
if (size >= 200)
numRepeats = 1;
else if (size >= 100)
numRepeats = 10;
if (size == 3) {
for (int k = 0; k < numRepeats; k++)
}
for (int k = 0; k < numRepeats; k++)
#ifdef BIAS_HAVE_OPENCV
if (useOpenCV)
else
#else
#endif
if (size <= MAX_PRINT_SIZE)
if (size == 3) {
}
if (size <= MAX_PRINT_SIZE)
cout << endl <<
"SVD S = " << svd.
GetS() << endl;
if (size == 3) {
cout << endl <<
"SVD3x3 S = " << svd3x3.
GetS() << endl
<< endl << "diff. S = " << dS << endl;
}
if (size <= MAX_PRINT_SIZE)
if (size == 3) {
}
for (int i = 0; i < size; i++)
cout << endl <<
"SVD |M - U*S*VT| = " << dM.
NormL2() << endl;
if (dM.
NormL2() > errorThreshold) {
BIASERR("SVD failed!");
}
if (size == 3) {
for (int i = 0; i < size; i++)
S[i][i] = svd3x3.
GetS()[i];
cout << endl <<
"SVD3x3 |M - U*S*VT| = " << dM.
NormL2() << endl;
if (dM.
NormL2() > errorThreshold) {
BIASERR("SVD3x3 failed!");
if (size <= MAX_PRINT_SIZE)
}
}
cout << endl <<
"SVD time = " << (0.001 * timer.
GetRealTime()) / numRepeats <<
" ms";
if (size == 3) {
cout << endl <<
"SVD3x3 time = " << (0.001 * timer3x3.
GetRealTime()) / numRepeats <<
" ms";
}
cout << endl << endl << "2.) Invert matrix with SVD" << endl;
#ifdef BIAS_HAVE_OPENCV
if (useOpenCV)
else
#else
#endif
dM = M * invM - I;
cout << endl <<
"SVD |M*inv(M) - I| = " << dM.
NormL2() << endl;
if (dM.
NormL2() > errorThreshold) {
BIASERR("Matrix inversion with SVD failed!");
if (size <= MAX_PRINT_SIZE)
}
if (size == 3) {
dM = M * invM - I;
cout << endl <<
"SVD3x3 |M*inv(M) - I| = " << dM.
NormL2() << endl;
if (dM.
NormL2() > errorThreshold) {
BIASERR("Matrix inversion with SVD3x3 failed!");
}
}
cout << endl << endl << "3.) Compute matrix square root with SVD" << endl;
if (size <= MAX_PRINT_SIZE)
if (size <= MAX_PRINT_SIZE)
dM = MTM - (sqrtMTM * sqrtMTM);
cout << endl <<
"SVD |MT*M - sqrt(MT*M)^2| = " << dM.
NormL2() << endl;
if (dM.
NormL2() > errorThreshold) {
BIASERR("Matrix square root failed!");
if (size <= MAX_PRINT_SIZE)
}
if (size == 3) {
sqrtMTM = svd3x3.
Sqrt(MTM);
dM = MTM - (sqrtMTM * sqrtMTM);
cout << endl <<
"SVD3x3 |MT*M - sqrt(MT*M)^2| = " << dM.
NormL2() << endl;
if (dM.
NormL2() > errorThreshold) {
BIASERR("Matrix square root failed!");
}
}
cout << endl;
return 0;
}