source: issm/branches/trunk-jpl-damage/src/c/toolkits/petsc/patches/MatInvert.cpp@ 12004

/*!\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);

        #if defined(_HAVE_MATLAB_) && defined(_SERIAL_)
                MatConvert(inv, MATSEQAIJ,MAT_REUSE_MATRIX,&inv);
        #else
                MatConvert(inv, MATMPIAIJ,MAT_REUSE_MATRIX,&inv);
        #endif
        
        /*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;

}