/*
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"

#include <mkl.h>

struct pardiso_matrix
{
	long n;
	long nz;
	bool symmetric;
	vector<int> ia; 
	vector<int> ja; 
	vector<double> a; 
};

// now from mkl.h
//#if defined(_WIN32) || defined(_WIN64)
//#define pardiso_ PARDISO
//#else
//#define PARDISO pardiso_
//#endif 
//extern "C" int PARDISO
//	(void *, int *, int *, int *, int *, int *,
//	double *, int *, int *, int *, int *, int *,
//	int *, double *, double *, int *);
//extern "C" int mkl_get_max_threads();

struct pardiso_struct
{
	void *pt[64];
	int maxfct;
	int mnum;
	int mtype;
	int n;
	double *a;
	int *ia;
	int *ja;
	int idum; //dummy not used by PARDISO when iparm(5-1) != 1
	int nrhs;
	int iparm[64];
	int msglvl;
	double ddum;
	int error;
	pardiso_struct()
	{
//		fill(pt, pt + 64, void(0)); does not work
		for (int i = 0; i < 64; ++i)
			pt[i] = 0;
		maxfct = 1;
		mnum = 1;
		nrhs = 1;
		fill(iparm, iparm + 64, 0);
		iparm[2] = mkl_get_max_threads();
	}
};

class ClassPARDISO : public LinearSolver
{	
	bool verbose;
public:
	ClassPARDISO();
	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)
	{
		pardiso_struct *id = (pardiso_struct*)L;
		int phase = 33;
		id->msglvl = 0; 
		PARDISO (id->pt, &(id->maxfct), &(id->mnum), &(id->mtype), &phase,
		&(id->n), id->a, id->ia, id->ja, &(id->idum), &(id->nrhs),
		id->iparm, &(id->msglvl), in, out, &(id->error));
		if (id->error != 0)
			throw SolverError::General(stringAndNumber("PARDISO solve error ", id->error));		
	}
	virtual void* copyMatrix(void *mat1);	
	virtual void sumMatrices(void *mat1, void* mat2, const double &alpha);	
	virtual bool isSymmetric(void *mat)
	{
		return ((pardiso_matrix*)mat)->symmetric;
	}
	virtual bool supportSymmetric() {return true;}
	
	virtual void clearMatrix(void *mat)
	{
		delete ((pardiso_matrix*)mat);
	}
	virtual void clearFactor(void *L)
	{
		pardiso_struct *id = (pardiso_struct*)L;
		int phase = -1;
		PARDISO (id->pt, &(id->maxfct), &(id->mnum), &(id->mtype), &phase,
		&(id->n), &(id->ddum), id->ia, id->ja, &(id->idum), &(id->nrhs),
		id->iparm, &(id->msglvl), &(id->ddum), &(id->ddum), &(id->error));
	 	delete id;	
	}	
	virtual void clear()
	{
	}
};

extern bool MUMPSpositive;
extern int MUMPSoutofcore;

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

void* ClassPARDISO::setMatrix(SparseMatrix &mat)
{
	pardiso_matrix *A = new pardiso_matrix;
	int dim = mat.size(); 
	if (mat.IsSymmetric)
		A->symmetric = true;
	else
		A->symmetric = false;
	A->n = dim;
	A->nz = mat.nnz;
	(A->ia).resize(dim + 1);
	(A->ja).resize(A->nz);
	(A->a).resize(A->nz);
	int pos = 0;
	for (int i = 0; i < dim; ++i)
	{
		A->ia[i] = pos + 1;
		for (map<int, double>::const_iterator j = mat[i].begin(); j != mat[i].end(); ++j)
		{
			A->ja[pos] = (*j).first + 1;
			A->a[pos] = (*j).second;	
			++pos;
		}
	}
	A->ia[dim] = mat.nnz + 1;
	mat.clearMemory();
	return A;
}

void ClassPARDISO::mulMatrixByVec(void *mat, double *in, double *out)
{
	pardiso_matrix &m = *(pardiso_matrix*)mat;	
	int n = m.n;
 	for (int i = 0; i < n; ++i)
		out[i] = 0.;
	for (int i = 0; i < n; i++) 
	{
		for (int jp = m.ia[i] - 1; jp < m.ia[i+1] - 1; jp++) 
		{
			int j   = m.ja[jp] - 1;
			double Aij = m.a[jp];	
			out[i] += Aij*in[j];
			if (m.symmetric && i != j)
				out[j] += Aij*in[i];
		}
	}
}

void* ClassPARDISO::factor(void *mat, bool ClearMatrix, const string &param)
{
	pardiso_struct *id = new pardiso_struct;
	pardiso_matrix &m = *(pardiso_matrix*)mat;
	if (m.symmetric)
	{
		if (MUMPSpositive)
			id->mtype = 2;
		else
			id->mtype = -2;
	}
	else	
  	id->mtype = 11;
//ordering metis (5), or pord (4), or AMD (0), AMF (2), QAMD (6)
	if (param == "metis")	
		id->iparm[1] = 2;
	else if (param == "md")
		id->iparm[6] = 0;
	else
		warning("no such ordering in PARDISO - ignored. Using METIS");
//streams  
	if (verbose)
	  id->msglvl = 1; 
	else
	  id->msglvl = 0; 
	if (MUMPSoutofcore == 1)
		id->iparm[59] = 2;
	else
		id->iparm[59] = 0;
	id->n = m.n;
	id->ia = &*m.ia.begin(); 
	id->ja = &*m.ja.begin();
	id->a = &*m.a.begin();

	int phase = 11;
	PARDISO (id->pt, &(id->maxfct), &(id->mnum), &(id->mtype), &phase,
	&(id->n), id->a, id->ia, id->ja, &(id->idum), &(id->nrhs),
	id->iparm, &(id->msglvl),  &(id->ddum),  &(id->ddum), &(id->error));
	if (id->error != 0)
		throw SolverError::General(stringAndNumber("PARDISO analysis error ", id->error));			
	phase = 22;
	PARDISO (id->pt, &(id->maxfct), &(id->mnum), &(id->mtype), &phase,
	&(id->n), id->a, id->ia, id->ja, &(id->idum), &(id->nrhs),
	id->iparm, &(id->msglvl),  &(id->ddum),  &(id->ddum), &(id->error));
	if (id->error != 0)
		throw SolverError::General(stringAndNumber("PARDISO factorization error ", id->error));			
// Can we clear matrix?	
	if (ClearMatrix) 
		clearMatrix(mat);
	return id;
}

void* ClassPARDISO::copyMatrix(void *mt1)
{
	pardiso_matrix &mat1 = *(pardiso_matrix*)mt1;
	pardiso_matrix *A = new pardiso_matrix;
// A = mat1	
	A->n = mat1.n;
	A->symmetric = mat1.symmetric;
	(A->ia).reserve(mat1.ia.size());
	(A->ja).reserve(mat1.ja.size());
	(A->a).reserve(mat1.a.size());
	(A->ia).assign(mat1.ia.begin(), mat1.ia.end());
	(A->ja).assign(mat1.ja.begin(), mat1.ja.end());
	(A->a).assign(mat1.a.begin(), mat1.a.end());
	return A;
}

void ClassPARDISO::sumMatrices(void *mt1, void* mt2, const double &alpha)
{
	pardiso_matrix *mat1 = (pardiso_matrix*)mt1;
	pardiso_matrix *mat2 = (pardiso_matrix*)mt2;
	if (mat1->symmetric != mat2->symmetric)
		throw SolverError::General("PARDISO:SumMatrices: matrices have different symmetry");
	//matK + s0 matC				
	//nnz of mat2 is subset of mat1
	//indeces are sorted
	for (int i = 0; i < mat1->n; ++i)
	{
		int start = mat1->ia[i] - 1;
		for (int ii = mat2->ia[i] - 1; ii < mat2->ia[i + 1] - 1; ++ii)
		{
			int j1 = -1;
			for (int k = start; k < mat1->ia[i + 1] - 1; ++k)
			{
				if (mat2->ja[ii] == mat1->ja[k])
				{
					j1 = k;
					break;
				}
				else if (mat2->ja[ii] < mat1->ja[k])
					break;
			}
		 	if (j1 == -1)
				throw SolverError::General("PARDISO:SumMatrices: no such nnz in matK");
			start = j1 + 1;
			mat1->a[j1] += mat2->a[ii]*alpha;
		}
	}
}

LinearSolver* make_ClassPARDISO()
{
	return new ClassPARDISO;
}