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


#undef __FUNCT__ 
#define __FUNCT__ "MatInvert"

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

#include "../../../shared/shared.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)throw ErrorException(__FUNCT__," trying to invert a non square matrix!");

	/*Create identitiy matrix: */
	identity=NewMat(M,N,&sparsity,NULL,NULL);
	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);
	MatMatSolve(factored_matrix,identity,inv);

	/*Assemble inverse: */
	MatAssemblyBegin(inv,MAT_FINAL_ASSEMBLY);
	MatAssemblyEnd(inv,MAT_FINAL_ASSEMBLY);

	#ifdef _SERIAL_
		MatConvert(inv, MATSEQAIJ,MAT_REUSE_MATRIX,&inv);
	#else
		MatConvert(inv, MATMPIAIJ,MAT_REUSE_MATRIX,&inv);
	#endif
	
	/*Free ressources:*/
	MatFree(&identity);
	VecFree(&diagonal);
	KSPDestroy(ksp);

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

}

