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

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

using namespace::std;

extern "C" {
int decomp_(int *ndim, int *n, double *a, double *cond, int *ipvt, double *work);
int solve_(int *ndim, int *n, double *a, double *b, int *ipvt);
}

int main(int argc, char* argv[])
{
	if (argc < 2)
	{
		printf("Usage: main dim\n");
		printf("Please specify matrix dimensions\n");
		return 1;
	}
	int dim;
	dim = atoi(argv[1]);
	Matrix A(dim, dim); 
	vector<double> B(dim), work(dim);
	vector<int> ipvt(dim);
	double cond;
	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++)
			B[i] = rand()/maxr;
	Matrix Acopy(A);
	vector<double> Bcopy(B);
	Timing factor;
	decomp_(&dim, &dim, &*A.begin(), &cond, &*ipvt.begin(), &*work.begin());
	double time = factor.time();
	cout << "time to factor A(" << dim << "," << dim << ") is " << time << " s" << endl;
	cout << "cond is " << cond << endl;
	if (cond + 1. == cond)
	{
		cout << "matrix is singular" << endl;
		return 1;
	}
	Timing back;
	solve_(&dim, &dim, &*A.begin(), &*B.begin(), &*ipvt.begin());
	time = back.time();
	cout << "time for back substitution is " << time << " s" << endl;
  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];
		}
    normsol += pow((Bcopy[i] - sum), 2);
  }
  cout << "residual is " << sqrt(normsol/normmat) << endl;	

	return 0;
}