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

The program reads the matrix and RHS in the Matrix Market format then run MUMPS.

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

#include <fstream>
#include <iomanip>
#include "matrices.h"

#include "dmumps_c.h"

// MUMPS does not have matrix vector product
void mulMatrixByVec(int n, int symmetric, vector<double> &a, vector<int> &irn, vector<int> &jcn, double *in, double *out)
{
 	for (int i = 0; i < n; ++i)
		out[i] = 0.;
	for (int k = 0; k < a.size(); k++) 
	{
		out[irn[k] - 1] += a[k]*in[jcn[k] - 1];
		if (symmetric && irn[k] != jcn[k])
			out[jcn[k] - 1] += a[k]*in[irn[k] - 1];
	}
}

int main(int argc, char **argv)
{	
	if (argc < 2)
	{
		cout << "run_mumps matrix rhs" << endl;
  	return 1;	
	}
	try
	{
// 1) Reading matrix
  SparseMatrix mat;
	ifstream in(argv[1]);
	if (!in)
	{
		cout << "cannot open file " << argv[1] << endl;
		return 2;
	}
	Timing t0;
 	mat.read(in);
  t0.write("Reading matrix");
	if (mat.NRows() != mat.NCols())
	{
		cout << "Matrix is not squared " << endl;
		return 3;
	}
// 2) converting matrix to the MUMPS representation and setting up MUMPS
// 2a) setting up MUMPS parameters and initialize it
	DMUMPS_STRUC_C id;
  id.job = -1; 
  id.par = 1; 
	if (mat.IsSymmetric)
// if you know that you matrix is positive definite, id.sym should be 1, then it is faster
//  id.sym = 1;
	  id.sym = 2;
	else
		id.sym = 0;
  id.comm_fortran = 1;
	Timing t1;
  dmumps_c(&id);
  t1.write("Initialization");
//streams    
  id.icntl[0] = 6; 
  id.icntl[1] = 6; 
  id.icntl[2] = 6; 
  id.icntl[3] = 3;
//ordering metis (5), or pord (4), or AMD (0), AMF (2), QAMD (6)	
  id.icntl[6] = 5;
  vector<int> irn;
  vector<int> jcn;
  vector<double> a;

// 2b) converting SparseMatrix to the MUMPS representation		
	Timing r1;
	id.n = mat.NRows();
	id.nz = mat.nnz;
	irn.resize(id.nz);
	jcn.resize(id.nz);
	a.resize(id.nz);
	int k = 0;
	for (int i = 0; i < id.n; ++i)
		for (map<int, double>::iterator j = mat[i].begin(); j != mat[i].end(); ++j)
		{
			irn[k] = i + 1;
			jcn[k] = (*j).first + 1;
			a[k] = (*j).second;
			++k;
		}
	id.irn = &*irn.begin(); 
	id.jcn = &*jcn.begin();
	id.a = &*a.begin(); 
	mat.clearMemory();
	r1.write("Converting matrix");
	Timing t2;
// 3) LU factorization	
  id.job = 1;
  dmumps_c(&id);
  t2.write("Analysis");
	Timing t3;
  id.job = 2;
  dmumps_c(&id);
  t3.write("Factorization");
	if (argc < 3)
		return 0;
// 4) Reading RHS and setting up MUMPS
	ifstream in2(argv[2]);
	if (!in2)
	{
		cout << "cannot open file " << argv[2] << endl;
		return 4;
	}
	Matrix rhs;
	rhs.read(in2);
	if (id.n != rhs.NRows())
	{
		cout << "RHS is not compatible with the matrix" << endl;
		return 4;
	}
	Matrix rhs_copy(rhs);
	id.rhs = &*rhs.begin();
	id.nrhs = rhs.NCols();
	id.lrhs = id.n;
	Timing t4;
// 5) Back substitution	
  id.job = 3;
  dmumps_c(&id);
  t4.write("Back substitution");

// 6) compute the norm and write the solution out
//solution in rhs (it is overwritten), the correct rhs is in rhs_copy	
	vector<double> out(id.n);
	for (int i = 0; i < rhs.NCols(); ++i)
	{
		mulMatrixByVec(id.n, id.sym, a, irn, jcn, rhs.column(i), &*out.begin());
		double sum = 0.;
		for (int j = 0; j < id.n; ++j)
			sum += pow(rhs_copy(j, i) - out[j], 2);
		cout << "norm Ax - b for RHS " << i + 1 << " is " << sqrt(sum) << endl;
	}
  cout << "Solution is in " << argv[2] << ".solution" << endl;
	ofstream outf((string(argv[2]) + ".solution").c_str());
	rhs.write(outf);

// 7) clean up
	Timing t5;	
  id.job = -2; 
  dmumps_c(&id);
  t5.write("Cleaning MUMPS");
	}
	catch (BaseError &t)
	{
		cout << "Error is encountered:" << endl;
		cout << t << endl;
		return 5;
	}
	catch (...)
	{
		cout << "Not known error is encountered." << endl;
		return 6;
	}
  return 0;	
}