/*
 *
 *  This file is part of MUMPS 4.8.0, built on Fri Jul 25 14:46:02 2008
 *
 */
/* $Id: c_example.c 5043 2008-07-18 08:56:02Z pcombes $ */
/* Example program using the C interface to the 
 * double real arithmetic version of MUMPS, dmumps_c.
 * We solve the system A x = RHS with
 *   A = diag(1 2) and RHS = [1 4]^T
 * Solution is [1 2]^T */

/*
 * Domen Stadler, 18.12.2011
 *
 * USAGE: mumps_solver matrixfile.mtx rhsfile.mtx solutionfile.mtx
 *
 * Source code was modified to read matrix (matrix market format) and
 * right hand side vector (matrix market format) and solve it. Currently 
 * it enables only one RHS vector at a time and symmetric matrices. 
 * Solution is written to the solution file.
 *
 * TODO: - parsing first lines from input files to determine if the matrix is not 
 *			symmetric and according to this, change the parameter id.sym
 *			(read MUMPS userguide for more instructions)
 *		 - implementation for several right hand side vector
 *			(currently I do not know how to do that exactly)
 *		 - implementation of other parameters for MUMPS
 *		 - test on other windows machines (I only tested it on win xp 32bit)
 */
#include <stdio.h>
#include "headers/mpi.h"
#include "headers/dmumps_c.h"
#define JOB_INIT -1
#define JOB_END -2
#define USE_COMM_WORLD -987654
#include <iostream>
#include <fstream>
#include <cstdlib>
#include <string>


#if defined(UNICODE)
	#define _tcout wcout
#else
	#define _tcout cout
#endif


int main(int argc, char* argv[])
{
	char* matFile = argv[1];
	char* rhsFile=argv[2];
	char* lhsFile = argv[3];
	int * iRn;
	int * jCn;
	double * mValues;
	int * iRHSn;
	int * jRHSn;
	double * rhsValues;
	bool matrixReaded=false;
	bool vectorReaded=false;
	int i,j, values ,nR,nC,nonZVmatrix, nLRHS ,n_RHS ,nz_RHS,rhs_old=0,rhs_temp=0 ;
	double value=0.0; 

	std::string line;
	std::ifstream matrixfile;
	std::ifstream vectorfile;
	std::ofstream solutionfile;

	//reading matrix
	std::cout<<"Reading matrix file "<<matFile<<std::endl;
	matrixfile.open(matFile);
	if(matrixfile.is_open())
	{
		getline(matrixfile,line);
		std::cout<<line<<std::endl;
		matrixfile>>nR>>nC>>nonZVmatrix;
		std::cout<<"Number of rows: "<<nR<<", number of columns: "<<nC<<", number of non zero Values: "<<nonZVmatrix;
		iRn = (int*) malloc(nonZVmatrix*sizeof(int));
		jCn = (int*) malloc(nonZVmatrix*sizeof(int));
		mValues =(double*) malloc(nonZVmatrix*sizeof(double));
		matrixReaded=true;
		for(i=0;i<nonZVmatrix;i++)
		{
			if(matrixfile.good())
			{
				matrixfile>>iRn[i]>>jCn[i]>>mValues[i];
			}
			else
			{
				std::cout<<"End of matrix file before it was filled"<<std::endl;
				free(iRn);
				free(jCn);
				free(mValues);
				free(matFile);
				free(rhsFile);
				free(lhsFile);
				matrixfile.close();
				return 11;
			}
		}
		matrixfile.close();
	}
	else
	{
		std::cout<<"File was not opened"<<std::endl;
		if(matrixReaded)
		{
			free(iRn);
			free(jCn);
			free(mValues);
		}
		free(matFile);
		free(rhsFile);
		free(lhsFile);
		return 10;
	}
	std::cout<<"Matrix readed"<<std::endl;
	//control of the matrix reading
	/*for(i=0;i<nonZVmatrix;i++)
	{
		std::cout<<"Row: "<<iRn[i]<<", Column: "<<jCn[i]<<", Value: "<<mValues[i]<<std::endl;
	}*/

	//reading RHS vector
	std::cout<<"Reading Right Hand Side vector file "<<rhsFile<<std::endl;
	vectorfile.open(rhsFile);
	if(vectorfile.is_open())
	{
		getline(vectorfile,line);
		std::cout<<line<<std::endl;
		vectorfile>>nLRHS>>n_RHS>>nz_RHS;
		std::cout<<"Length of vector: "<<nLRHS<<", number of vectors: "<<n_RHS<<", number of non zero Values: "<<nz_RHS<<std::endl;
		rhsValues =(double*) malloc(nLRHS*sizeof(double));
		vectorReaded=true;
		for(i=0;i<nLRHS;i++)
		{rhsValues[i]=0.0;}
		std::cout<<"here"<<std::endl;
		double value;
		for(i=0;i<nz_RHS;i++)
		{
			if(vectorfile.good())
			{
				vectorfile>>j>>rhs_temp>>value;
				rhsValues[j-1]=value;
			}
			else
			{
				std::cout<<"End of vector file before it was filled"<<std::endl;
				free(iRn);
				free(jCn);
				free(mValues);
				free(iRHSn);
				free(jRHSn);
				free(rhsValues);
				vectorfile.close();
				return 21;
			}
		}
		vectorfile.close();
	}
	else
	{
		std::cout<<"File was not opened"<<std::endl;
		if(matrixReaded)
		{
			free(iRn);
			free(jCn);
			free(mValues);
		}
		if(vectorReaded)
		{
			free(iRHSn);
			free(jRHSn);
			free(rhsValues);
		}
		free(matFile);
		free(rhsFile);
		free(lhsFile);
		return 20;
	}
	std::cout<<"Right Hand Side vector readed"<<std::endl;
	//control of the vector reading
	/*for(i=0;i<nLRHS;i++)
	{
		std::cout<<"Row: "<<i<<", Vector: "<<0<<", Value: "<<rhsValues[i]<<std::endl;
	}*/


	DMUMPS_STRUC_C id;

	int myid, ierr;
#if defined(MAIN_COMP)
	argv = &name;
#endif
	std::cout<<"Setting MPI environment"<<std::endl;
	ierr = MPI_Init(&argc, &argv);
	ierr = MPI_Comm_rank(MPI_COMM_WORLD, &myid);

	/* Initialize a MUMPS instance. Use MPI_COMM_WORLD */
	std::cout<<"Initializing MUMPS instance"<<std::endl;
	id.job=JOB_INIT; id.par=1; id.sym=2;id.comm_fortran=USE_COMM_WORLD;
	dmumps_c(&id);
	/* Define the problem on the host */
#define ICNTL(I) icntl[(I)-1] /* macro s.t. indices match documentation */
	std::cout<<"Filling MUMPS arrays"<<std::endl;
	if (myid == 0) {
		id.n = nR; id.nz =nonZVmatrix; id.irn=iRn; id.jcn=jCn;
		id.a = mValues; id.rhs = rhsValues;
	}
	else
	{
		std::cout<<"Could not fill the MUMPS arrays"<<std::endl;
		if(matrixReaded)
		{
			free(iRn);
			free(jCn);
			free(mValues);
		}
		if(vectorReaded)
		{
			free(iRHSn);
			free(jRHSn);
			free(rhsValues);
		}
		free(matFile);
		free(rhsFile);
		free(lhsFile);
		return 30;
	}
	/* No outputs */
	id.ICNTL(1)=-1; id.ICNTL(2)=-1; id.ICNTL(3)=-1; id.ICNTL(4)=0;
	/* Call the MUMPS package. */
	id.job=6;
	std::cout<<"Solving the linear system"<<std::endl;
	dmumps_c(&id);
	std::cout<<"Terminating the MUMPS instance"<<std::endl;
	id.job=JOB_END; dmumps_c(&id); /* Terminate instance */
	if (myid == 0) 
	{
		std::cout<<"Writing solution to the file"<<lhsFile<<std::endl;
		solutionfile.open(lhsFile);
		solutionfile<<"%%MatrixMarket matrix coordinate real general" <<std::endl;
		solutionfile<<nLRHS<<" 1 "<<nLRHS <<std::endl;
		for(i=0;i<nLRHS;i++)
		{
			solutionfile<<" "<<(i+1)<<" 1 "<<rhsValues[i]<<std::endl;
		}
		solutionfile.close();
		std::cout<<"Solution written"<<std::endl;
	}
	else
	{
		std::cout<<"Could not write the solution file"<<std::endl;
		if(matrixReaded)
		{
			free(iRn);
			free(jCn);
			free(mValues);
		}
		if(vectorReaded)
		{
			free(iRHSn);
			free(jRHSn);
			free(rhsValues);
		}
		free(matFile);
		free(rhsFile);
		free(lhsFile);
		return 31;
	}
	std::cout<<"Finalizing MPI"<<std::endl;
	ierr = MPI_Finalize();

	std::cout<<"Freeing allocated arrays"<<std::endl;
	if(matrixReaded)
	{
		free(iRn);
		free(jCn);
		free(mValues);
	}
	if(vectorReaded)
	{
		free(iRHSn);
		free(jRHSn);
		free(rhsValues);
	}
	free(matFile);
	free(rhsFile);
	free(lhsFile);
	return 0;
}