/*
Copyright (C) 2006-2009 Evgenii Rudnyi, http://Evgenii.Rudnyi.Ru/

The program reads the matrix and RHS in the Matrix Market format then run a sparse solver to solve A x = rhs.

This software is a copyrighted work licensed under the terms, described
in the file "FREE_LICENSE". 
*/

#include <fstream>
#include "solvers.h"

string solver = "MUMPS";
string parameter = "metis";
bool verboseSolver = false;
string matrix;
string rhsFile;

extern bool MUMPSpositive;
extern int MUMPSoutofcore;

const char* HELP ="\n\
    run_sparse_solver [-options] matrix [rhs]\n\
    -v  print information from solvers \n\
    -s  solver (MUMPS, TAUCSlltmf, TAUCSllt, TAUCSldlt, TAUCSlltll, TAUCSlltooc, TAUCSlu, UMFPACK PARDISO) \n\
    -p  parameter (TAUCS: metis, genmmd, md, mmd, amd) \n\
                  (MUMPS: metis, AMD, AMF, QAMD) \n\
                  (MUMPS,PARDISO: out-of-core, ldlt) \n";

void readArguments(int argc, char *argv[])
{
	string par;
	int ii;
  for (int i = 1; i < argc; i++)
  {
    if (argv[i][0] == '-')
    {
      switch (argv[i][1])
      {
        case 'v' :
					verboseSolver = true;
          break;
        case 's' :
          par = argv[i] + 2;
          if (par.empty() && i + 1 < argc && argv[i + 1][0] != '-')
            par = argv[++i];
					if (!par.empty()) solver = par;
					cout << "INPUT: Solver is " << solver << "." << endl;
          break;
       	case 'p' :
          par = argv[i] + 2;
          if (par.empty() && i + 1 < argc && argv[i + 1][0] != '-')
            par = argv[++i];
					if (!par.empty()) 
					{	
						if (par == "ldlt")
						{
							MUMPSpositive = false;
							cout << "INPUT: Forcing LDLT for MUMPS and PARDISO."<< endl;
						}
						else if (par == "out-of-core")
						{
							MUMPSoutofcore = 1;
							cout << "INPUT: Forcing out-of-core for MUMPS and PARDISO."<< endl;
						}
						else
						{
							parameter = par;
							cout << "INPUT: Parameter is " << parameter << "." << endl;
						}
					}
          break;
       	default:
          cout << "WARNING: Wrong option " << argv[i] << ": ignored." << endl; 
          break;
      }
    }
    else
    {
			 ifstream test(argv[i]);
			 if (!test)
				 cout << "WARNING: File " << argv[i] << " does not exist: ignored." << endl;
			 else
			 {
				 if (matrix.empty())
					 matrix = argv[i];
				 else if (rhsFile.empty())
					 rhsFile = argv[i];
				 else
				   cout << "WARNING: I expect only two arguments. Argument " << argv[i] << " is ignored." << endl;
			 }
    }
	}
}


int main(int argc, char *argv[])
{
	if (argc < 2)
	{
		cout << HELP << endl;
		return 1;
	}
	try
	{
		readArguments(argc, argv);
		if (matrix.empty())
		{
			cout << "Nothing to do." << endl << endl;				
			cout << HELP << endl;
			return 2;
		}
// 1) Initializing solver
		LinearSolver *sol = LinearSolver::create(solver); 
		if (verboseSolver)
			sol->setVerbose();

// 2) Reading matrix
		SparseMatrix mat;
		ifstream in(matrix.c_str());
		Timing t2;
		mat.read(in);
		t2.write("Reading matrix");
		if (mat.NRows() != mat.NCols())
		{
			cout << "Matrix is not squared " << endl;
			return 3;
		}
		int dim = mat.NRows();

// 3) converting matrix to the solver representation
    Timing t3;
    if (mat.IsSymmetric && !(sol->supportSymmetric()))
      mat.makeNotSymmetric();			
    void *m = sol->setMatrix(mat);
		t3.write("Converting matrix");

// 4) LU factorization
    Timing t4;
    void *L = sol->factor(m, false, parameter);
    t4.write("Factorization");
  	if (rhsFile.empty())
	  	return 4;

// 5) Reading RHS
		ifstream in2(argv[2]);
		Matrix rhs;
		rhs.read(in2);
		if (dim != rhs.NRows())
		{
			cout << "RHS is not compatible with the matrix" << endl;
			return 4;
		}

// 6) Back substitution	
    Timing t6;
		Matrix x(rhs.NRows(), rhs.NCols());
		for (int i = 0; i < rhs.NCols(); ++i)
			sol->mulInverseByVec(L, rhs.column(i), x.column(i));
    t6.write("Back substitution");		

// 7) compute the norm and write the solution out
		vector<double> out(dim);
		for (int i = 0; i < rhs.NCols(); ++i)
		{
			sol->mulMatrixByVec(m, x.column(i), &*out.begin());
			double sum = 0.;
			for (int j = 0; j < dim; ++j)
				sum += pow(rhs(j, i) - out[j], 2);
			cout << "norm Ax - b for RHS " << i + 1 << " is " << sqrt(sum) << endl;
		}
		cout << "Solution is in " << rhsFile << ".solution" << endl;
		ofstream outf((rhsFile + ".solution").c_str());
		x.write(outf);  

// 8) clean up
    sol->clearFactor(L);
		sol->clearMatrix(m);
		delete sol;
	}
	catch (BaseError &t)
	{
		cout << "Error is encountered during solution:" << endl;
		cout << t << endl;
		return 5;
	}
	catch (...)
	{
		cout << "ERROR: Not known error during Arnoldi." << endl;
		return 5;
	}
	return 0;
}