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

class TAUCS : public LinearSolver
{
public:
	virtual void* setMatrix(SparseMatrix &mat);
	virtual void setVerbose()
	{
		taucs_logfile("stdout");
	}
	virtual void mulMatrixByVec(void *mat, double *in, double *out)
	{
		taucs_dccs_times_vec((taucs_ccs_matrix*)mat, in, out);
	}

	virtual void* copyMatrix(void *mat1);	
	virtual void sumMatrices(void *mat1, void* mat2, const double &alpha);	
	virtual bool isSymmetric(void *mat)
	{
		return ((taucs_ccs_matrix*)mat)->flags & TAUCS_SYMMETRIC;
	}
	
	virtual void clearMatrix(void *mat)
	{
		taucs_dccs_free((taucs_ccs_matrix*)mat);
	}
	
};

void* TAUCS::setMatrix(SparseMatrix &mat)
{
	taucs_ccs_matrix *A;
	int dim = mat.size(); 
	A = taucs_dccs_create(dim, dim, mat.nnz);
	A->flags = TAUCS_DOUBLE;
	if (mat.IsSymmetric)
		A->flags |= TAUCS_SYMMETRIC | TAUCS_LOWER;
	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);				
	}
	int pos = 0;
	for (int i = 0; i < dim; ++i)
	{
		A->colptr[i] = pos;
		for (map<int, double>::const_iterator j = mat[i].begin(); j != mat[i].end(); ++j)
		{
			A->rowind[pos] = (*j).first;
			A->values.d[pos] = (*j).second;	
			++pos;
		}
	}
	A->colptr[dim] = mat.nnz;
	mat.clearMemory();
	return A;
}

void* TAUCS::copyMatrix(void *mt1)
{
	taucs_ccs_matrix *A;
	taucs_ccs_matrix *mat1 = (taucs_ccs_matrix*)mt1;
	int dim = mat1->n; 
// A = mat1	
	A = taucs_dccs_create(dim, dim, mat1->colptr[dim]);
	copy(mat1->colptr, mat1->colptr + dim + 1, A->colptr);
	copy(mat1->rowind, mat1->rowind + mat1->colptr[dim], A->rowind);
	copy(mat1->values.d, mat1->values.d + mat1->colptr[dim], A->values.d);
	A->flags = mat1->flags;
	return A;
}	

void TAUCS::sumMatrices(void *mt1, void* mt2, const double &alpha)
{
	taucs_ccs_matrix *mat1 = (taucs_ccs_matrix*)mt1;
	taucs_ccs_matrix *mat2 = (taucs_ccs_matrix*)mt2;
//matK + s0 matC	
	if ((mat1->flags & TAUCS_SYMMETRIC)  && !(mat2->flags & TAUCS_SYMMETRIC))
		throw SolverError::General("SumMatrices: matrix K is symmetric and matrix C is not");
//to the end of the loop	
			if ((mat2->flags & TAUCS_SYMMETRIC) && 
					!(mat1->flags & TAUCS_SYMMETRIC) //&& 
//					(i != mat2->rowind[ii])
					)
			{
				throw SolverError::General("SumMatrices: matC symm and matK not symm: not implemented yet");
			}						
//nnz of mat2 is subset of mat1
//indeces are sorted
	int dim = mat1->n; 
	for (int i = 0; i < dim; ++i)
	{
		int start = mat1->colptr[i];
		for (int ii = mat2->colptr[i]; ii < mat2->colptr[i + 1]; ++ii)
		{
			int j1 = -1;
			for (int k = start; k < mat1->colptr[i + 1]; ++k)
			{
				if (mat2->rowind[ii] == mat1->rowind[k])
				{
					j1 = k;
					break;
				}
				else if (mat2->rowind[ii] < mat1->rowind[k])
					break;
			}
		 	if (j1 == -1)
				throw SolverError::General("SumMatrices: TAUCS: no such nnz in matK");
			start = j1 + 1;
			mat1->values.d[j1] += mat2->values.d[ii]*alpha;
		}
	}
}

class TAUCSdirectLLT : public TAUCS
{	
	int* perm;
	int* invperm;
	vector<double> vec_tmp;
	int dim;

	virtual void* factoring(taucs_ccs_matrix* mat) = 0;
	virtual void substitution(void* L, double *sol, double *x) = 0;

public:
	TAUCSdirectLLT() {perm = 0;}
	virtual void* factor(void *mat, bool ClearMatrix, const string &param);
	virtual void mulInverseByVec(void *L, double *in, double *out)
	{
		taucs_dvec_permute(dim, in, &*vec_tmp.begin(), perm);
		substitution(L, out, &*vec_tmp.begin());
  	taucs_dvec_ipermute(dim, out, &*vec_tmp.begin(), perm);
		for (int i = 0; i < dim; ++i)
			out[i] = vec_tmp[i];
	}
	virtual bool supportSymmetric() {return true;}
	virtual void clear()
	{
		if (perm)
		{
			free(perm);
			perm = 0;
		}
	}
};

void* TAUCSdirectLLT::factor(void *mt, bool ClearMatrix, const string &param)
{
	taucs_ccs_matrix *mat = (taucs_ccs_matrix*)mt;
	dim = mat->n;
	if (!(mat->flags & TAUCS_SYMMETRIC))
		throw SolverError::TAUCS("Matrix to factor is not symmetric: try another solver");
	
	void *L;
	taucs_ccs_matrix* ccs_temp;
	Timing reod;
	taucs_ccs_order(mat, &perm, &invperm, (char*)param.c_str());
	ccs_temp = taucs_dccs_permute_symmetrically(mat, perm, invperm);
	if (ClearMatrix) 
		taucs_dccs_free((taucs_ccs_matrix*)mat);
	free(invperm);
	reod.write("Reodering");
	
	L = factoring(ccs_temp);
	taucs_dccs_free(ccs_temp);
	if (!L)
		throw SolverError::TAUCS("Factoring is impossible: try another solver");
	vec_tmp.resize(dim);
	return L;
}

class TAUCSllt : public TAUCSdirectLLT
{
	virtual void* factoring(taucs_ccs_matrix* mat)
	{
		return taucs_dccs_factor_llt(mat, 0., 0);
	}
	virtual void substitution(void* L, double *sol, double *x)
	{
		taucs_dccs_solve_llt(L, sol, x);
	}
public:	
	virtual void clearFactor(void *L)
	{
		taucs_dccs_free((taucs_ccs_matrix*)L);
	}
};

LinearSolver* make_TAUCSllt()
{
	return new TAUCSllt;
}

class TAUCSldlt : public TAUCSdirectLLT
{
	virtual void* factoring(taucs_ccs_matrix* mat)
	{
		return taucs_dccs_factor_ldlt(mat);
	}
	virtual void substitution(void* L, double *sol, double *x)
	{
		taucs_dccs_solve_ldlt(L, sol, x);
	}
public:	
	virtual void clearFactor(void *L)
	{
		taucs_dccs_free((taucs_ccs_matrix*)L);
	}
};

LinearSolver* make_TAUCSldlt()
{
	return new TAUCSldlt;
}

class TAUCSllt_mf : public TAUCSdirectLLT
{
	virtual void* factoring(taucs_ccs_matrix* mat)
	{
		return taucs_dccs_factor_llt_mf(mat);
	}
	virtual void substitution(void* L, double *sol, double *x)
	{
		taucs_dsupernodal_solve_llt(L, sol, x);
	}
public:	
	virtual void clearFactor(void *L)
	{
		taucs_dsupernodal_factor_free(L);
	}
};

LinearSolver* make_TAUCSllt_mf()
{
	return new TAUCSllt_mf;
}

class TAUCSllt_ll : public TAUCSdirectLLT
{
	virtual void* factoring(taucs_ccs_matrix* mat)
	{
		return taucs_dccs_factor_llt_ll(mat);
	}
	virtual void substitution(void* L, double *sol, double *x)
	{
		taucs_dsupernodal_solve_llt(L, sol, x);
	}
public:	
	virtual void clearFactor(void *L)
	{
		taucs_dsupernodal_factor_free(L);
	}
};

LinearSolver* make_TAUCSllt_ll()
{
	return new TAUCSllt_ll;
}

class TAUCSllt_ooc : public TAUCSdirectLLT
{
	virtual void* factoring(taucs_ccs_matrix* mat)
	{
		void* L = taucs_io_create_multifile("tmp_llt_ooc");
		double mem = taucs_available_memory_size();
		information(stringAndNumber("INFO: memory requested is ", mem));
		int ret = taucs_dooc_factor_llt(mat, (taucs_io_handle*)L, mem);
		if (ret == -1) 
			throw SolverError::TAUCS("Factoring is impossible: try another solver");
		return L;
	}
	virtual void substitution(void* L, double *sol, double *x)
	{
		taucs_dooc_solve_llt(L, sol, x);
	}
public:	
	virtual void clearFactor(void *L)
	{
		taucs_io_delete((taucs_io_handle*)L);
	}	
};

LinearSolver* make_TAUCSllt_ooc()
{
	return new TAUCSllt_ooc;
}

class TAUCSlu : public TAUCS
{	
	taucs_io_handle *L;
	int* perm;
	int* invperm;
public:	

	virtual void* factor(void *mt, bool clearMatrix, const string &param)
	{
		taucs_ccs_matrix *mat = (taucs_ccs_matrix*)mt;
		if (mat->flags & TAUCS_SYMMETRIC)
			throw SolverError::TAUCS("LU: Matrix to factor is symmetric: use symmetric solver");
		Timing reod;
		taucs_ccs_order(mat, &perm, &invperm, "colamd");
		free(invperm);
		reod.write("Reodering");

		L = taucs_io_create_multifile("tmp_lu_ooc");
		double mem = taucs_available_memory_size();
		information(stringAndNumber("INFO: memory requested is ", mem));
		int ret = taucs_ooc_factor_lu(mat, perm, L, mem);
		if (clearMatrix) taucs_dccs_free(mat);
		free(perm);
		if (ret == -1) 
			throw SolverError::TAUCS("Factoring is impossible: try another solver");
		return L;
	}
	virtual void mulInverseByVec(void *L, double *in, double *out)
	{
		taucs_dooc_solve_lu((taucs_io_handle*)L, out, in);
	}
	virtual bool supportSymmetric() {return false;}
	virtual void clearFactor(void *L)
	{
		taucs_io_delete((taucs_io_handle*)L);
	}	
};

LinearSolver* make_TAUCSlu()
{
	return new TAUCSlu;
}