// Evgenii B. Rudnyi, http://Evgenii.Rudnyi.Ru

#include <iostream>
#include <stdlib.h>
#include <math.h>
#include "matrix.h"

using namespace::std;

int main(int argc, char* argv[])
{
	if (argc < 2)
	{
		printf("Usage: main dim nrs\n");
		printf("Please specify matrix dimensions\n");
		return 1;
	}
	int dim, nrhs;
	dim = atoi(argv[1]);
	if (argc < 3)
		nrhs = 1;
	else
		nrhs = atoi(argv[2]);
	Matrix A(dim, dim); 
	Matrix B(dim, nrhs);;
	vector<int> ipvt(dim);
	srand(86456);
	double maxr=(double)RAND_MAX;
	for (int i = 0; i < dim; i++)
		for (int j = 0; j < dim; j++)
			A(i, j) = rand()/maxr;
	for (int i = 0; i < dim; i++)
		for (int j = 0; j < nrhs; j++)
			B(i, j) = rand()/maxr;
	Matrix Acopy(A);
	Matrix Bcopy(B);
	Timing factor;
	int info;
	info = dgetf2(A, ipvt);
	double time = factor.time();
	cout << "time to factor A(" << dim << "," << dim << ") with dgetf2 is " << time << " s" << endl;
	cout << "info is " << info << endl;
	A = Acopy;
	Timing factor2;
	info = dgetrf(A, ipvt);
	time = factor2.time();
	cout << "time to factor A(" << dim << "," << dim << ") with dgetrf is " << time << " s" << endl;
	cout << "info is " << info << endl;
	if (info != 0)
	{
		cout << "something went wrong" << endl;
		return 1;
	}
	Timing back;
	info = dgetrs(A, B, ipvt);
	time = back.time();
	cout << "time for back substitution is " << time << " s" << endl;
	cout << "info is " << info << endl;
	if (info != 0)
	{
		cout << "something went wrong" << endl;
		return 1;
	}
  double normmat = 0.;
	double normsol = 0.;
  for (int i = 0; i < dim; ++i)
  {
    double sum = 0.;
    for (int j = 0; j < dim; ++j)
		{
			normmat += Acopy(i, j)*Acopy(i, j);
      sum += Acopy(i, j)*B(j, 0);
		}
    normsol += pow((Bcopy(i, 0) - sum), 2);
  }
	for (int k = 1; k < nrhs; ++k)
		for (int i = 0; i < dim; ++i)
  	{
    	double sum = 0.;
	    for (int j = 0; j < dim; ++j)
	      sum += Acopy(i, j)*B(j, k);
    	normsol += pow((Bcopy(i, k) - sum), 2);
  	}
  cout << "residual is " << sqrt(normsol/normmat) << endl;	
	return 0;
}