/*!\file:  MatInvert.cpp
 * \brief  invert petsc matrix using LU factorization, and multiple right hand sides.
 */ 

/*Petsc includes: */
#include "petscmat.h"
#include "petscvec.h"
#include "petscksp.h"

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

void MatInvert(Mat* pinv, Mat matrix){


	/*output: */
	Mat inv=NULL;

	/*intermediary: */
	int M,N;
	double sparsity=.001;

	Mat factored_matrix=NULL;
	Mat identity=NULL;
	Vec diagonal=NULL;
	KSP ksp=NULL;
	PC pc=NULL;

	/*Some checks: */
	MatGetSize(matrix,&M,&N);
	if(M!=N) _error_("trying to invert a non square matrix!");

	/*Create identitiy matrix: */
	identity=NewMat(M,N,sparsity);
	diagonal=NewVec(M);
	VecSet(diagonal,1.0);
	MatDiagonalSet(identity,diagonal,INSERT_VALUES);
	MatAssemblyBegin(identity,MAT_FINAL_ASSEMBLY);
	MatAssemblyEnd(identity,MAT_FINAL_ASSEMBLY);
	MatConvert(identity, MATMPIDENSE,MAT_REUSE_MATRIX,&identity);
	
	/*Initialize inverse: */
	MatDuplicate(identity,MAT_DO_NOT_COPY_VALUES,&inv);

	/* Compute X=inv(A) by MatMatSolve(), from tests in petsc test27 */
	KSPCreate(PETSC_COMM_WORLD,&ksp);
	KSPSetOperators(ksp,matrix,matrix,DIFFERENT_NONZERO_PATTERN);
	KSPGetPC(ksp,&pc);
	PCSetType(pc,PCLU);
	KSPSetUp(ksp);
	//PCGetFactoredMatrix(pc,&factored_matrix); //not found any replacement yet
	MatMatSolve(factored_matrix,identity,inv);

	/*Assemble inverse: */
	MatAssemblyBegin(inv,MAT_FINAL_ASSEMBLY);
	MatAssemblyEnd(inv,MAT_FINAL_ASSEMBLY);
		
	MatConvert(inv, MATMPIAIJ,MAT_REUSE_MATRIX,&inv);
	
	/*Free ressources:*/
	MatFree(&identity);
	VecFree(&diagonal);
	#if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
	KSPDestroy(ksp);
	#else
	KSPDestroy(&ksp);
	#endif

	/*Assign output pointers:*/
	*pinv=inv;

}

