/*
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 <umfpack.h>
}

struct umfpack_matrix
{
	long n;
	long nz;
	vector<UF_long> Ap; //colptr
	vector<UF_long> Ai; //rowind
	vector<double> Ax;  //value.d
};

class UMFPACK : public LinearSolver
{	
	vector<double> control;
	vector<double> info;
	bool verbose;
public:
	UMFPACK() : control(UMFPACK_CONTROL), info(UMFPACK_INFO) 
	{
		umfpack_dl_defaults(&*control.begin());
		control[UMFPACK_IRSTEP] = 0;
		control[UMFPACK_PRL] = 2;
		verbose = false;
	}
	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)
	{
		int o = umfpack_dl_solve(UMFPACK_A, 0, 0, 0, out, in, L, &*control.begin(), &*info.begin());
//		if (verbose)
//			umfpack_dl_report_info(&*control.begin(), &*info.begin());
		if (o != 0)
			umfpack_dl_report_status(&*control.begin(), o);
	}
	virtual void* copyMatrix(void *mat1);	
	virtual void sumMatrices(void *mat1, void* mat2, const double &alpha);	
	virtual bool isSymmetric(void *mat)
	{
		return false;
	}
	virtual bool supportSymmetric() {return false;}
	
	virtual void clearMatrix(void *mat)
	{
		delete ((umfpack_matrix*)mat);
	}
	virtual void clearFactor(void *L)
	{
		umfpack_dl_free_numeric(&L);
	}	
	virtual void clear()
	{
	}
};

void* UMFPACK::setMatrix(SparseMatrix &mat)
{
	umfpack_matrix *A = new umfpack_matrix;
	int dim = mat.size(); 
	if (mat.IsSymmetric)
		throw SolverError::UMFPACK("Matrix to transform is symmetric");
	else
	{
		SparseMatrix tmp(dim);
		for (int i = 0; i < dim; ++i)
			for (map<int, double>::const_iterator j = mat[i].begin(); j != mat[i].end(); ++j)
				tmp[(*j).first][i] = (*j).second;
		tmp.nnz = mat.nnz;
		mat.swap(tmp);				
	}
	A->n = dim;
	A->nz = mat.nnz;
	(A->Ap).resize(dim + 1);
	(A->Ai).resize(A->nz);
	(A->Ax).resize(A->nz);
	int pos = 0;
	for (int i = 0; i < dim; ++i)
	{
		A->Ap[i] = pos;
		for (map<int, double>::const_iterator j = mat[i].begin(); j != mat[i].end(); ++j)
		{
			A->Ai[pos] = (*j).first;
			A->Ax[pos] = (*j).second;	
			++pos;
		}
	}
	A->Ap[dim] = mat.nnz;
	mat.clearMemory();
	return A;
}

void UMFPACK::mulMatrixByVec(void *mat, double *in, double *out)
{
	umfpack_matrix *m = (umfpack_matrix*)mat;	
	long n = m->n;
 	for (long i = 0; i < n; ++i)
		out[i] = 0.;
	for (long j = 0; j < n; j++) 
	{
		for (long ip = (m->Ap)[j]; ip < (m->Ap)[j+1]; ip++) 
		{
			long i   = (m->Ai)[ip];
			double Aij = (m->Ax)[ip];	
			out[i] += Aij*in[j];
		}
	}
}

void* UMFPACK::factor(void *mat, bool ClearMatrix, const string &param)
{
	void *L, *S;
	umfpack_matrix *m = (umfpack_matrix*)mat;	
	int out = umfpack_dl_symbolic(m->n, m->n, &*(m->Ap.begin()), &*(m->Ai.begin()), &*(m->Ax.begin()), &S,  &*control.begin(), &*info.begin());
//		if (verbose)
//			umfpack_di_report_info(&*control.begin(), &*info.begin());
	if (out != 0)
	{
		umfpack_dl_report_status(&*control.begin(), out);
		throw SolverError::UMFPACK(stringAndNumber("UMFPACK: symbolic analysis is not possible. Code is ", out));
	}
	double mem = 0;
	mem = atof(param.c_str());
//	cout << "info[UMFPACK_SYMMETRIC_DMAX] " << info[UMFPACK_SYMMETRIC_DMAX] << endl;
	if (mem != 0)
	{
		information(string("UMFPACK: UMFPACK_ALLOC_INIT is now ").append(ObjToString(-mem)));
		control[UMFPACK_ALLOC_INIT] = -mem;
		control[UMFPACK_FRONT_ALLOC_INIT] = -mem;
	}
	out = umfpack_dl_numeric(&*(m->Ap.begin()), &*(m->Ai.begin()), &*(m->Ax.begin()), S, &L,  &*control.begin(), &*info.begin());
	if (verbose)
		umfpack_dl_report_info(&*control.begin(), &*info.begin());
	if (out != 0)
	{
		umfpack_dl_report_status(&*control.begin(), out);
		throw SolverError::UMFPACK(stringAndNumber("UMFPACK: numeric analysis is not possible. Code is ", out));
	}
	umfpack_dl_free_symbolic(&S);
	if (ClearMatrix) 
		clearMatrix(mat);
	return L;
}

void* UMFPACK::copyMatrix(void *mt1)
{
	umfpack_matrix *A = new umfpack_matrix;
	umfpack_matrix *mat1 = (umfpack_matrix*)mt1;
	int dim = mat1->n; 
// A = mat1	
	A->n = dim;
	A->nz = mat1->Ap[dim];
	(A->Ap).reserve(dim + 1);
	(A->Ai).reserve(A->nz);
	(A->Ax).reserve(A->nz);
	(A->Ap).assign(mat1->Ap.begin(), mat1->Ap.end());
	(A->Ai).assign(mat1->Ai.begin(), mat1->Ai.end());
	(A->Ax).assign(mat1->Ax.begin(), mat1->Ax.end());
	return A;
}

void UMFPACK::sumMatrices(void *mt1, void* mt2, const double &alpha)
{
	umfpack_matrix *mat1 = (umfpack_matrix*)mt1;
	umfpack_matrix *mat2 = (umfpack_matrix*)mt2;
//matK + s0 matC				
//nnz of mat2 is subset of mat1
//indeces are sorted
	int dim = mat1->n; 
	for (int i = 0; i < dim; ++i)
	{
		int start = mat1->Ap[i];
		for (int ii = mat2->Ap[i]; ii < mat2->Ap[i + 1]; ++ii)
		{
			int j1 = -1;
			for (int k = start; k < mat1->Ap[i + 1]; ++k)
			{
				if (mat2->Ai[ii] == mat1->Ai[k])
				{
					j1 = k;
					break;
				}
				else if (mat2->Ai[ii] < mat1->Ai[k])
					break;
			}
		 	if (j1 == -1)
				throw SolverError::General("SumMatrices: no such nnz in matK");
			start = j1 + 1;
			mat1->Ax[j1] += mat2->Ax[ii]*alpha;
		}
	}
}

LinearSolver* make_UMFPACK()
{
	return new UMFPACK;
}