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

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

#include "solvers.h"

extern "C" {
#include "dmumps_c.h"
}

struct mumps_matrix
{
	bool symmetric;
	int n;
// used to be DMUMPS_INT in 4.7.3	and now MUMPS_INT in 4.8.3
  vector<MUMPS_INT> irn;
  vector<MUMPS_INT> jcn;
  vector<double> a;
};

class MUMPS : public LinearSolver
{	
	bool verbose;
public:
	MUMPS();
	virtual void* setMatrix(SparseMatrix &mat);
	virtual void setVerbose()
	{
		verbose = true;
	}
	virtual void mulMatrixByVec(void *mat, double *in, double *out);

	virtual void* factor(void *mat, bool ClearMatrix, const string &param);
	virtual void mulInverseByVec(void *L, double *in, double *out)
	{
		DMUMPS_STRUC_C *id = (DMUMPS_STRUC_C*)L;
		for (int i = 0; i < id->n; ++i)
			out[i] = in[i];
		id->rhs = out;
	  id->job = 3;
		//streams    
	  id->icntl[0] = 0; 
  	id->icntl[1] = 0; 
  	id->icntl[2] = 0; 
  	id->icntl[3] = 0;
  	dmumps_c(id);
		if (id->infog[0] > 0)
			warning(stringAndNumber("MUMPS solve warning ", id->infog[1]));
		else if (id->infog[0] < 0)
			throw SolverError::General(stringAndNumber("MUMPS solve error ", id->infog[1]));		
	}
	virtual void* copyMatrix(void *mat1);	
	virtual void sumMatrices(void *mat1, void* mat2, const double &alpha);	
	virtual bool isSymmetric(void *mat)
	{
		return ((mumps_matrix*)mat)->symmetric;
	}
	virtual bool supportSymmetric() {return true;}
	
	virtual void clearMatrix(void *mat)
	{
		delete ((mumps_matrix*)mat);
	}
	virtual void clearFactor(void *L)
	{
		DMUMPS_STRUC_C *id = (DMUMPS_STRUC_C*)L;
// JOB_END		
	  id->job = -2; 
  	dmumps_c(id);
	 	delete id;	
	}	
	virtual void clear()
	{
	}
};

bool MUMPSpositive = true;
int MUMPSoutofcore = 0;

MUMPS::MUMPS()
{
	verbose = false;
}

void* MUMPS::setMatrix(SparseMatrix &mat)
{
	mumps_matrix *A = new mumps_matrix;
	int dim = mat.size(); 
	if (mat.IsSymmetric)
		A->symmetric = true;
	else
		A->symmetric = false;
	A->n = dim;
	(A->irn).resize(mat.nnz);
	(A->jcn).resize(mat.nnz);
	(A->a).resize(mat.nnz);
	int k = 0;
	for (int i = 0; i < A->n; ++i)
		for (map<int, double>::iterator j = mat[i].begin(); j != mat[i].end(); ++j)
		{
			(A->irn)[k] = i + 1;
			(A->jcn)[k] = (*j).first + 1;
			(A->a)[k] = (*j).second;
			++k;
		}	
	mat.clearMemory();
	return A;
}

void MUMPS::mulMatrixByVec(void *mat, double *in, double *out)
{
	mumps_matrix &m = *(mumps_matrix*)mat;	
	int n = m.n;
 	for (int i = 0; i < n; ++i)
		out[i] = 0.;
	for (int k = 0; k < m.a.size(); k++) 
	{
		out[m.irn[k] - 1] += m.a[k]*in[m.jcn[k] - 1];
		if (m.symmetric && m.irn[k] != m.jcn[k])
			out[m.jcn[k] - 1] += m.a[k]*in[m.irn[k] - 1];
	}
}

void* MUMPS::factor(void *mat, bool ClearMatrix, const string &param)
{
	DMUMPS_STRUC_C *id = new DMUMPS_STRUC_C;
	mumps_matrix &m = *(mumps_matrix*)mat;
  id->par=1; 
  id->comm_fortran=1;
	if (m.symmetric)
	{
		if (MUMPSpositive)
			id->sym = 1;
		else
			id->sym = 2;
	}
	else	
  	id->sym = 0;
	// JOB_INIT; 
  id->job=-1;
  dmumps_c(id);
//ordering metis (5), or pord (4), or AMD (0), AMF (2), QAMD (6)
	if (param == "metis")	
		id->icntl[6] = 5;
	else if (param == "pord")
		id->icntl[6] = 4;
	else if (param == "AMD")
		id->icntl[6] = 0;
	else if (param == "AMF")
		id->icntl[6] = 2;
	else if (param == "QAMD")
		id->icntl[6] = 6;
	else
	{
		id->icntl[6] = 5;
		warning("no such ordering in MUMPS - ignored. Using METIS");
	}
//streams  
	if (verbose)
	{	
	  id->icntl[0] = 6; 
  	id->icntl[1] = 6; 
  	id->icntl[2] = 6; 
  	id->icntl[3] = 3;
	}
	else
	{	
	  id->icntl[0] = 0; 
  	id->icntl[1] = 0; 
  	id->icntl[2] = 0; 
  	id->icntl[3] = 0;
	}
	id->icntl[21] = MUMPSoutofcore;
	if (MUMPSoutofcore == 1)
	{
//		for (int i = 0; i < 150; ++i)
//			id->ooc_tmpdir[i] = 0;
		id->ooc_tmpdir[0] = '.';
		id->ooc_tmpdir[1] = 0;
	}
	id->n = m.n;
  id->nz = m.a.size();
	id->irn = &*m.irn.begin(); 
	id->jcn = &*m.jcn.begin();
	id->a = &*m.a.begin();

  id->job = 1;
  dmumps_c(id);
	if (id->infog[0] > 0)
		warning(stringAndNumber("MUMPS analysis warning ", id->infog[1]));
	else if (id->infog[0] < 0)
		throw SolverError::General(stringAndNumber("MUMPS analysis error ", id->infog[1]));			
  id->job = 2;
  dmumps_c(id);
	if (id->infog[0] > 0)
		warning(stringAndNumber("MUMPS factoring warning ", id->infog[1]));
	else if (id->infog[0] < 0)
		throw SolverError::General(stringAndNumber("MUMPS factoring error ", id->infog[1]));			
// Can we clear matrix?	
	if (ClearMatrix) 
		clearMatrix(mat);
	return id;
}

void* MUMPS::copyMatrix(void *mt1)
{
	mumps_matrix &mat1 = *(mumps_matrix*)mt1;
	mumps_matrix *A = new mumps_matrix;
// A = mat1	
	A->n = mat1.n;
	A->symmetric = mat1.symmetric;
	(A->irn).reserve(mat1.irn.size());
	(A->jcn).reserve(mat1.jcn.size());
	(A->a).reserve(mat1.a.size());
	(A->irn).assign(mat1.irn.begin(), mat1.irn.end());
	(A->jcn).assign(mat1.jcn.begin(), mat1.jcn.end());
	(A->a).assign(mat1.a.begin(), mat1.a.end());
	return A;
}

void MUMPS::sumMatrices(void *mt1, void* mt2, const double &alpha)
{
	mumps_matrix &mat1 = *(mumps_matrix*)mt1;
	mumps_matrix &mat2 = *(mumps_matrix*)mt2;
	if (mat1.symmetric != mat2.symmetric)
		throw SolverError::General("MUMPS:SumMatrices: matrices have different symmetry");
	//matK + s0 matC				
	int start = 0;
	for (int i = 0; i < mat2.a.size(); ++i)
	{
//cout << "mat2 " << mat2.irn[i] << " " << mat2.jcn[i] << endl;	
		int j1 = -1;
//it could be that we are still on the previous row in mat1
		for (int k = start; mat1.irn[k] < mat2.irn[i]; ++k)
			++start;
		for (int k = start; mat1.irn[k] == mat2.irn[i]; ++k)
		{
//cout << " mat1 " << mat1.irn[k] << " " << mat1.jcn[k] << endl;	
			if (mat1.jcn[k] == mat2.jcn[i])
			{
				j1 = k;
				break;
			}
		}
		if (j1 == -1)
			throw SolverError::General("SumMatrices: no such nnz in matK");
		start = j1 + 1;
		mat1.a[j1] += mat2.a[i]*alpha;
	}
}

LinearSolver* make_MUMPS()
{
	return new MUMPS;
}