/*!\file SolverxPetsc
 * \brief Petsc implementation of solver
 */

#include "./Solverx.h"
#include "../../shared/shared.h"
#include "../../include/include.h"
#include "../../io/io.h"

#ifdef HAVE_CONFIG_H
	#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

void	SolverxPetsc(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters){

	/*Output: */
	Vec        uf               = NULL;

	/*Intermediary: */
	int        local_m,local_n,global_m,global_n;
	int        analysis_type;

	/*Solver */
	KSP        ksp              = NULL;
	PC         pc               = NULL;
	int        iteration_number;
	int        solver_type;
	bool       fromlocalsize    = true;
	#if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
	PetscTruth flag,flg;
	#else
	PetscBool flag,flg;
	#endif

	/*Stokes: */
	IS         isv=NULL;
	IS         isp=NULL;

	#if _PETSC_MAJOR_ >= 3 
	char ksp_type[50];
	#endif


	/*Display message*/
	_printf_(VerboseModule(),"   Solving\n");
	#if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
	if(VerboseSolver())PetscOptionsPrint(stdout);
	#else
	if(VerboseSolver())PetscOptionsView(PETSC_VIEWER_STDOUT_WORLD);
	#endif

	/*First, check that f-set is not NULL, i.e. model is fully constrained: {{{*/
	_assert_(Kff);
	MatGetSize(Kff,&global_m,&global_n); _assert_(global_m==global_m);
	if(!global_n){
		*puf=NULL; return;
	}
	/*}}}*/
	/*Initial guess logic here: {{{1*/
	/*Now, check that we are not giving an initial guess to the solver, if we are running a direct solver: */
	#if _PETSC_MAJOR_ >= 3 
	PetscOptionsGetString(PETSC_NULL,"-ksp_type",ksp_type,49,&flg);
	if (strcmp(ksp_type,"preonly")==0)uf0=NULL;
	#endif

	/*If initial guess for the solution exists, use it to create uf, otherwise, 
	 * duplicate the right hand side so that the solution vector has the same structure*/
	if(uf0){
		VecDuplicate(uf0,&uf); VecCopy(uf0,uf);
	}
	else{
		MatGetLocalSize(Kff,&local_m,&local_n);uf=NewVec(local_n,fromlocalsize);
	}
	/*}}}*/
	/*Process petsc options to see if we are using special types of external solvers: {{{1*/
	PetscOptionsDetermineSolverType(&solver_type);

	/*In serial mode, the matrices have been loaded as MPIAIJ or AIJ matrices. 
	 We need to convert them if we are going to run the solvers successfully: */
	#ifdef _SERIAL_
	#if _PETSC_MAJOR_ == 2 
	if (solver_type==MUMPSPACKAGE_LU){
		/*Convert Kff to MATTAIJMUMPS: */
		MatConvert(Kff,MATAIJMUMPS,MAT_REUSE_MATRIX,&Kff);
	}
	if (solver_type==MUMPSPACKAGE_CHOL){
		/*Convert Kff to MATTSBAIJMUMPS: */
		MatConvert(Kff,MATSBAIJMUMPS,MAT_REUSE_MATRIX,&Kff);
	}
	if (solver_type==SPOOLESPACKAGE_LU){
		/*Convert Kff to MATTSBAIJMUMPS: */
		MatConvert(Kff,MATAIJSPOOLES,MAT_REUSE_MATRIX,&Kff);
	}
	if (solver_type==SPOOLESPACKAGE_CHOL){
		/*Convert Kff to MATTSBAIJMUMPS: */
		MatConvert(Kff,MATSBAIJSPOOLES,MAT_REUSE_MATRIX,&Kff);
	}
	if (solver_type==SUPERLUDISTPACKAGE){
		/*Convert Kff to MATTSBAIJMUMPS: */
		MatConvert(Kff,MATSUPERLU_DIST,MAT_REUSE_MATRIX,&Kff);
	}
	if (solver_type==StokesSolverEnum){
		_error_("Petsc 2 does not support multi-physics solvers");
	}
	#endif
	#endif
	/*}}}*/
	/*Check the solver is available: {{{1*/
	if(solver_type==MUMPSPACKAGE_LU || solver_type==MUMPSPACKAGE_CHOL){
	#if _PETSC_MAJOR_ >=3
		#ifndef _HAVE_MUMPS_
		_error_("requested MUMPS solver, which was not compiled into ISSM!\n");
		#endif

	#endif
	}
	/*}}}*/
	/*Prepare solver:{{{1*/
	KSPCreate(MPI_COMM_WORLD,&ksp);
	KSPSetOperators(ksp,Kff,Kff,DIFFERENT_NONZERO_PATTERN);
	KSPSetFromOptions(ksp);

	#if defined(_SERIAL_) && _PETSC_MAJOR_==3
	/*Specific solver?: */
	KSPGetPC(ksp,&pc);
	if (solver_type==MUMPSPACKAGE_LU){
		#if _PETSC_MINOR_==1
		PCFactorSetMatSolverPackage(pc,MAT_SOLVER_MUMPS);
		#else
		PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS);
		#endif
	}
	#endif

	#if defined(_PARALLEL_) && _PETSC_MAJOR_==3
	/*Stokes: */
	if (solver_type==StokesSolverEnum){
		/*Make indices out of doftypes: */
		if(!df)_error_("need doftypes for Stokes solver!\n");
		DofTypesToIndexSet(&isv,&isp,df,StokesSolverEnum);

		/*Set field splits: */
		KSPGetPC(ksp,&pc);
		#if _PETSC_MINOR_==1
		PCFieldSplitSetIS(pc,isv);
		PCFieldSplitSetIS(pc,isp);
		#else
		PCFieldSplitSetIS(pc,PETSC_NULL,isv);
		PCFieldSplitSetIS(pc,PETSC_NULL,isp);
		#endif

	}
	#endif

	/*}}}*/
	/*If there is an initial guess for the solution, use it, except if we are using the MUMPS direct solver, where any initial guess will crash Petsc: {{{1*/
	if (uf0){
		if( (solver_type!=MUMPSPACKAGE_LU) && (solver_type!=MUMPSPACKAGE_CHOL) && (solver_type!=SPOOLESPACKAGE_LU) && (solver_type!=SPOOLESPACKAGE_CHOL) && (solver_type!=SUPERLUDISTPACKAGE)){
			KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
		}
	}
	/*}}}*/
	
	if(VerboseSolver())KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);

	/*Solve: */
	KSPSolve(ksp,pf,uf);
	
	/*Check convergence*/
	KSPGetIterationNumber(ksp,&iteration_number);
	if (iteration_number<0) _error_("%s%i"," Solver diverged at iteration number: ",-iteration_number);

	/*Free resources:*/
	KSPFree(&ksp);
		
	/*Assign output pointers:*/
	*puf=uf;
}