[petsc-users] Matrix inverse and LU decomposition

[Bartłomiej W]

That solved the problem, thank you very much!
Afterall, it ocurred that I don't need this procedure since I will only use diagonal matrices and inverse of such are trivial.

It will come in handy someday though, I am sure.
Thanks again.

--- 30.3.11 (Śr), Hong Zhang <hzhang at mcs.anl.gov> napisał(a):

Od: Hong Zhang <hzhang at mcs.anl.gov>
Temat: Re: [petsc-users] Matrix inverse and LU decomposition
Do: "PETSc users list" <petsc-users at mcs.anl.gov>
Data: 30 Marzec 2011 (Środa), 15:07

You need call
MatSetType(I,MATDENSE);MatSetType(Dinv,MATDENSE);
See petsc-dev/src/mat/examples/tests/ex1.c (attached) on how to useMatMatSolve().Can you send me your complete code so I can check why Dinv is empty?

Hong

On Wed, Mar 30, 2011 at 8:31 AM, Bartłomiej W <bartlomiej.wach at yahoo.pl> wrote:

Hello, 

I want to obtain the inverse of a matrix, even though it's not advised by the FAQ ;-) and I have a problem that earlier responses didn't solve for me. Simply speaking : the inverse matrix that is a result of  DI=Dinv as told by the FAQ, is empty for me. Here is my code(the significant part) :


  ierr = MatCreate(PETSC_COMM_WORLD,&D);CHKERRQ(ierr);
  ierr = MatSetSizes(D,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
  ierr = MatSetFromOptions(D);CHKERRQ(ierr);

  ierr = MatCreate(PETSC_COMM_WORLD,&Dinv);CHKERRQ(ierr);

  ierr = MatSetSizes(Dinv,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
  ierr = MatSetFromOptions(Dinv);CHKERRQ(ierr);

  ierr = MatCreate(PETSC_COMM_WORLD,&I);CHKERRQ(ierr);
  ierr =
 MatSetSizes(I,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
  ierr = MatSetFromOptions(I);CHKERRQ(ierr);

(...)

  ierr = MatAssemblyEnd(D,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(I,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

  ierr = MatAssemblyEnd(Dinv,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

  MatView(D,PETSC_VIEWER_STDOUT_SELF);
    
  ierr = MatGetOrdering(D,  MATORDERING_NATURAL,  &perm,  &iperm); CHKERRQ(ierr);     
  ierr = MatFactorInfoInitialize(&info); CHKERRQ(ierr);

  ierr = MatLUFactor(D, perm, iperm, &info); CHKERRQ(ierr);

  MatView(D,PETSC_VIEWER_STDOUT_SELF);
  MatView(I,PETSC_VIEWER_STDOUT_SELF);
    
  ierr = MatMatSolve(D,I,Dinv); CHKERRQ(ierr);

  MatView(Dinv,PETSC_VIEWER_STDOUT_SELF);


// Results in such output:

//====
 MATRIX D
row 0: (0, 2) 
row 1: (1, 1) 
row 2: (2, 1) 
row 3: (3, 1) 

//==== MATRIX D after LU decomposition
row 0: (0, 0.5) 
row 1: (1, 1) 
row 2: (2, 1) 
row 3: (3, 1) 

//==== Identity matrix

row 0: (0, 1) 
row 1: (1, 1) 
row 2: (2, 1) 
row 3: (3, 1) 

//==== inverse of D  (result of MatMatSolve(D,I,Dinv))
row 0:
row 1:
row 2:
row 3:


The approach with :

    ierr = MatGetOrdering(D,  MATORDERING_NATURAL,  &perm,  &iperm); CHKERRQ(ierr);     

    ierr = MatGetFactor(D, MAT_SOLVER_PETSC, MAT_FACTOR_LU,&F);
    ierr = MatFactorInfoInitialize(&info); CHKERRQ(ierr);
    ierr = MatLUFactorSymbolic(F,D,perm,iperm,&info);CHKERRQ(ierr);
    ierr = MatLUFactorNumeric(F,D,&info);CHKERRQ(ierr);

    ierr = MatLUFactor(D,
 perm, iperm, &info); CHKERRQ(ierr);

Provides same output

Thank You for Your help.



-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20110330/f00750a0/attachment.htm>