#include "stdafx.h"
#include <fstream>
#include <iostream>
#include "smatrix.h"
#include "RndSmatrix.h"
#include "LUDecomp.h"
#include "fastheap.h"
#include <iostream>
#include <time.h>
using namespace sm;
smatrix<double> CoefInit5x5()
{
smatrix<double> c;
c(1, 2) = 1;
c(2, 1) = -2; c(2, 3) = 7; c(2, 5) = 1;
c(3, 3) = -9;
c(4, 1) = -6; c(4, 4) = 6; c(4, 5) = 4;
c(5, 1) = 7; c(5, 2) = -6; c(5, 4) = 6;
return(c);
}
double My5x5[] = {
0, 1, 0, 0, 0,
-2, 0, 7, 0, 1,
0, 0, -9, 0, 0,
-6, 0, 0, 6, 4,
7, -6, 0, 6, 0
};
smatrix<double> ArrayInit5x5()
{
return(smatrix<double>(My5x5, 5, 5));
}
smatrix<double> InsertInit5x5()
{
smatrix<double> c(5, 5, 11);
double *d = My5x5;
for(size_t i = 1; i <= 5; i++)
{
for(size_t j = 1; j <= 5; j++, d++)
{
if(*d != 0)
c.insert(i, j, *d);
}
}
return(c);
}
template <class t> smatrix<t> MyMultiply(const smatrix<t> &a,
const smatrix<t> &b)
{
smatrix<t> c;
for(size_t i = 1; i <= a.iDim; i++)
{
for(size_t j = 1; j <= b.jDim; j++)
c(i, j) = a.row(i) * b.col(j);
}
return(c);
}
template <class t> smatrix<t> MyAdd(const smatrix<t> &a,
const smatrix<t> &b)
{
smatrix<t> c;
for(size_t i = 1; i <= a.iDim; i++)
c.row(i) = a.row(i) + b.row(i);
return(c);
}
template <class t> smatrix<t> MyUpperTrans(const smatrix<t> &m)
{
smatrix<t> c;
for(size_t i = 1; i <= m.iDim; i++)
c(rng(i, m.iDim), i) = m(i, rng(i, m.iDim));
return(c);
}
size_t TestLUDecomp(ostream &os, size_t Size,
size_t NZ, size_t xNZ, size_t Iterations)
{
os << "Test Linear Solutions Using LU Decomposition ";
os << endl;
EnableFastHeap();
srand(0);
time_t Total(0), Start, Stop;
for(size_t n = 0; n < Iterations; n++)
{
// Construct a random linear system, Ax = b
TRndSolvable<double, double> A(Size, -10, 10, NZ);
os << "A" << endl << A << endl;
smatrix<double> x;
if(xNZ == Size)
x = TRndDMatrix<double, double>(Size, 1, -10, 10);
else
x = TRndSMatrix<double, double>(Size, 1, -10,
10, xNZ);
os << "x" << endl << x << endl;
smatrix<double> b = A * x;
os << "b" << endl << b << endl;
// Solve for 'x' given 'A' and 'b'
time(&Start);
TLUDecomp<double> LU(A, true);
smatrix<double> xCalc = LU.Solve(b);
time(&Stop);
Total += Stop - Start;
// Print how close the 'xCalc' is to 'x'
os << "Residual = " << (xCalc - x).maxCoef << endl;
os << endl << endl;
}
DisableFastHeap();
return(Total / Iterations);
}
size_t TestLUInverse(ostream &os, size_t Size,
size_t NZ, size_t xNZ, size_t Iterations)
{
os << "Testing Inverse Using LU Decomposition " << endl;
EnableFastHeap();
srand(0);
time_t Total(0), Start, Stop;
for(size_t n = 0; n < Iterations; n++)
{
smatrix<double> A = TRndDMatrix<double, double>(
Size, Size, -10, 10);
time(&Start);
TLUDecomp<double> LU(A, true);
smatrix<double> InvA = LU.Inverse();
time(&Stop);
Total += Stop - Start;
smatrix<double> Ident = smIdentity<double>(Size);
os << InvA << endl;
os << "Residual = " << (Ident - (A * InvA)).maxCoef;
os << endl << endl << endl;
}
DisableFastHeap();
return(Total / Iterations);
}
template <class t> smatrix<t> JonesMayerInv(
const smatrix<t> &m, size_t Iterations)
{
smatrix<t> a(m);
t Scaler = t(1.0) / m.FrobeniousNormal() * 1.1;
a *= Scaler;
smatrix<t> x(a);
x.Transpose();
for(size_t n = 0; n < Iterations; n++)
x = (x * 2) - ((x * a) * x);
x *= Scaler;
return(x);
}
void TestJonesMayerInv(ostream &os, size_t Size,
size_t Iterations)
{
os << "Testing Jones and Mayer Inversion" << endl;
EnableFastHeap();
srand(0);
for(size_t n = 0; n < Iterations; n++)
{
smatrix<double> A = TRndDMatrix<double, double>(
Size, Size, -10, 10);
smatrix<double> InvA = JonesMayerInv(A, 25);
smatrix<double> Ident = smIdentity<double>(Size);
os << InvA << endl;
os << "Residual = " << (Ident - (A * InvA)).maxCoef;
os << endl << endl << endl;
}
DisableFastHeap();
}
void PrintSolvability(ostream &os, smatrix<double> &m)
{
TLUDecomp<double> LU(m, true);
os << m << endl;
os << (!LU ? "Unsolvable" : "Solvable");
os << endl << endl;
}
void Unsolvable(ostream &os)
{
double _m[] = {
1.0, 1.0 - .00000001,
1.0 + .00000001, 1.0
};
smatrix<double> m(_m, 2, 2);
PrintSolvability(os, m);
}
void Solvable(ostream &os)
{
double _m[] = {
1.0, 1.0 - .0000001,
1.0 + .0000001, 1.0
};
smatrix<double> m(_m, 2, 2);
PrintSolvability(os, m);
}
void TestDeterminant(ostream &os, size_t Size,
size_t Iterations)
{
os << "Testing Determinant Calculation" << endl;
EnableFastHeap();
srand(0);
for(size_t n = 0; n < Iterations; n++)
{
double lnScaledDet;
smatrix<double> A = TRndDMatrix<double, double>(
Size, Size, -10, 10);
smatrix<double> temp(A);
int Sign = temp.lnScaledDet(lnScaledDet, true);
double Scaler = exp(lnScaledDet / Size);
smatrix<double> Unscaled = A * (1.0 / Scaler);
if(Sign < 0)
os << "Residual = " << (1.0 + Unscaled.det(true));
else
os << "Residual = " << (1.0 - Unscaled.det(true));
os << endl;
}
DisableFastHeap();
}
size_t smatrix<double>::PrintWidth = 8;
size_t smatrix<double>::PrintPrecision = 3;
int main(int argc, char* argv[])
{
cout << ArrayInit5x5() << endl;
cout << CoefInit5x5() << endl;
cout << InsertInit5x5() << endl;
fstream Fig1("Figure1.txt", ios::out);
smatrix<double> a = ArrayInit5x5();
a.PrintListing(Fig1);
cout << (a * 2.0) << endl;
cout << (2.0 * a) << endl;
smatrix<double> b = MyUpperTrans(a);
cout << b << endl;
cout << (a + b) << endl;
cout << MyAdd(a, b) << endl;
cout << (a * b) << endl;
cout << MyMultiply(a, b) << endl;
a = TRndDMatrix<double, double>(4, 5, -6, 7);
cout << a << endl;
b = TRndSMatrix<double, int>(8, 9, -11, 12, 50);
cout << b << endl;
TestJonesMayerInv(cout, 10, 10);
TestLUDecomp(cout, 10, 50, 10, 10);
TestLUInverse(cout, 10, 50, 10, 10);
Solvable(cout);
Unsolvable(cout);
cout << endl;
TestDeterminant(cout, 10, 50);
cout << "Average time: " <<
TestLUDecomp(cout, 2000, 15000, 2000, 10)
<< " seconds" << endl;
return 0;
}