[petsc-users] Inaccurate Matrix Inversion using MatMatSolve

[Barry Smith]

  Please send the code. This should work fine.

  Barry


> On Jun 25, 2020, at 2:43 PM, Rakesh Halder <rhalder at umich.edu> wrote:
> 
> Hi all,
> 
> I'm using PETSc for matrix calculations as part of a model order reduction code. An algorithm I'm using requires that I compute the explicit inverse of a matrix, as it needs to be used in matrix-matrix products. 
> 
> The matrix is small (I'll show a 5x5 example), and I first find the LU factorization of it as follows:
> 
>   MatGetOrdering(A,MATORDERINGNATURAL,&perm,&iperm);
>   MatFactorInfoInitialize(&info);
>   MatLUFactor(A,perm,iperm,&info);
> 
> I then find the inverse of it:
> 
>   MatMatSolve(A,I,B);
> 
> Where I is the identity matrix, and B is the inverse. The subroutine does calculate an inverse, although not accurately; the matrix A I looked at is:
> 
> 5.23E-02	1.86E-02	2.67E-02	4.58E-02	3.55E-02
> 6.37E-03	5.86E-02	5.07E-03	1.64E-02	1.36E-02
> -4.07E-02	3.99E-03	5.50E-02	1.77E-02	-3.21E-02
> 1.96E-02	-5.53E-03	4.02E-02	-5.37E-02	1.80E-02
> 1.51E-02	1.70E-02	-1.57E-02	-3.75E-03	-5.64E-02
> 
> And when I take the product A*B, I don't recover the identity matrix, but rather:
> 
> 9.97E-01	2.77E-04	-1.66E-03	1.25E-04	-1.10E-03
> -6.79E-04	9.99E-01	-1.72E-04	6.24E-04	1.16E-03
> 2.70E-03	-3.98E-04	9.98E-01	-1.51E-04	6.82E-04
> -3.35E-04	-4.82E-04	-1.03E-03	9.99E-01	9.80E-05
> -9.50E-04	-8.70E-04	8.37E-04	-3.68E-04	9.99E-01
> 
> I believe this is causing large inaccuracies in my program, as the diagonal and off-diagonal entries have large errors associated with them. I am wondering if there is a way to perhaps tighten the tolerance of the subroutine, or if there are other methods I can use. I would like to avoid using KSP and solving for the inverse vector-by-vector if possible.
> 
> Thanks,
> 
> Rakesh Halder

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200625/db752f3b/attachment-0001.html>