/*!\file Solverx
 * \brief solver
 */

#include "./Solverx.h"

#undef __FUNCT__ 
#define __FUNCT__ "Solverx"

#include "../shared/shared.h"

void	Solverx( Vec* puf, Mat Kff, Vec pf, Vec uf0, char* solver_string){

	/*output: */
	Vec uf=NULL;

	/*intermediary: */
	int local_m,local_n;
	
	/*Solver*/
	KSP ksp=NULL;
	PC  pc=NULL;
	int iteration_number;
	PetscTruth flag;
	int solver_type;

	/*If initial guess for solution exists, use it to create uf, otherwise, 
	 * duplicate right hand side so that solution vector has same structure*/
	if(uf0){
		VecDuplicate(uf0,&uf); VecCopy(uf0,uf);
	}
	else{
		MatGetLocalSize(Kff,&local_m,&local_n);uf=NewVecFromLocalSize(local_n);
	}

	/*Before preparing the solver, add options to the options database*/
	PetscOptionsInsertMultipleString(solver_string);

	/*Process solver_string to see if we are not using special types of external solvers: */
	PetscOptionsDetermineSolverType(&solver_type,solver_string);

	/*obsolete? : */
//	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);
//	}
//

	/*Prepare solver*/
	KSPCreate(MPI_COMM_WORLD,&ksp);
	KSPSetOperators(ksp,Kff,Kff,DIFFERENT_NONZERO_PATTERN);
	KSPSetFromOptions(ksp);

	/*specific solver?: */
	KSPGetPC(ksp,&pc);
	if (solver_type==MUMPSPACKAGE_LU){
		PCFactorSetMatSolverPackage(pc,MAT_SOLVER_MUMPS);
	}

	/*If initial guess for solution, use it, except if we are using the MUMPS direct solver, where any initial 
	 * guess will crash Petsc: */
	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);
		}
	}
		

	#ifdef _ISSM_DEBUG_ 
		KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);
		KSPGetInitialGuessNonzero(ksp,&flag); 
		_printf_("User provided initial guess solution: %i\n",flag);
	#endif

	KSPSolve(ksp,pf,uf);
	
	/*Check convergence*/
	KSPGetIterationNumber(ksp,&iteration_number);
	if (iteration_number<0){
		throw ErrorException(__FUNCT__,exprintf("%s%i"," Solver diverged at iteration number: ",-iteration_number));
	}
	else{
		#ifdef _ISSM_DEBUG_ 
		_printf_("%s%i%s\n","Solver converged after ",iteration_number," iterations");
		#endif
	}

	/*Free ressources:*/
	KSPFree(&ksp);

	
	/*Assign output pointers:*/
	*puf=uf;
}