Cholesky factorization problem

[Hong Zhang]

In petsc, we do A = R*D*R^T, in which D
acturally stores inv(d_i).
When A is in blocked form, i.e., bs=2 in your case,
d_i is a bsxbs block,
D holds inverse(d_i) which is computed by LAPACK.
Our MatView() only displays R without diagonal blocks
because this is the data structure used by MatSolve().

Why do you want to see R? For varifying the correctness
of the code or its numerical values?

Hong

On Fri, 13 Apr 2007 JUNWANG at uwm.edu wrote:

> Hi, Hong:
>
>     Really appreciate your help!
>     The MatSolve() works fine. It provides the right solution. How about if i
> wanna see the elements in the factorized up-triangle matix R, where A =R'R.
> Since the MatView() cannot show the logical matrix elements, which function
> should i call?
>     Thank you !
>
> Jun
>
> Quoting Hong Zhang <hzhang at mcs.anl.gov>:
>
> >
> > Jun,
> >
> > Your code looks fine. The factored matrix
> > "result" is an internal matrix to be used by MatSolve,
> > thus, its actural storage doesn't match the logical
> > matrix elements. MatView() only display off-diagonal blocks
> > of the factored matrix (these values may not match the logical
> > factored matrix values!).
> >
> > We'll add a note to MatView() when user display
> > an internal factored matrix.
> >
> > If you get wrong MatSolve() with your code,
> > let us know.
> >
> > Thanks for reporting if,
> >
> > Hong
> >
> >
> > On Fri, 13 Apr 2007 junwang at uwm.edu wrote:
> >
> > >
> > >
> > > Hello:
> > >    I have written a test sample program for using the Cholesky
> > Factorization
> > > function, here is the code i wrote, it compiled successfully. But the
> > result is
> > > not correct, is anyone who can help me out? Thank you!
> > > ***************************************************************************
> > >
> > > #include "petsc.h"
> > > #include "petscmat.h"
> > > #include <iostream.h>
> > > #include <stdlib.h>
> > >
> > > int main(int argc,char **args)
> > > {
> > >   PetscErrorCode ierr;
> > >   Mat            mat;
> > >   int ind1[2],ind2[2];
> > >   PetscScalar temp[2*2];
> > >   PetscInt *nnz=new PetscInt[3];
> > >
> > >   PetscInitialize(&argc,&args,0,0);
> > >   nnz[0]=2;nnz[1]=1;nnz[1]=1;
> > >
> > >   ierr = MatCreateSeqSBAIJ(PETSC_COMM_SELF,2,6,6,0,nnz,&mat);CHKERRQ(ierr);
> > >   ind1[0]=0;ind1[1]=1;
> > >   temp[0]=3;temp[1]=2;temp[2]=0;temp[3]=3;
> > >   ierr = MatSetValues(mat,2,ind1,2,ind1,temp,INSERT_VALUES);CHKERRQ(ierr);
> > >   ind2[0]=4;ind2[1]=5;
> > >   temp[0]=1;temp[1]=1;temp[2]=2;temp[3]=1;
> > >   ierr = MatSetValues(mat,2,ind1,2,ind2,temp,INSERT_VALUES);CHKERRQ(ierr);
> > >   ind1[0]=2;ind1[1]=3;
> > >   temp[0]=4;temp[1]=1;temp[2]=1;temp[3]=5;
> > >   ierr = MatSetValues(mat,2,ind1,2,ind1,temp,INSERT_VALUES);CHKERRQ(ierr);
> > >   ind1[0]=4;ind1[1]=5;
> > >   temp[0]=5;temp[1]=1;temp[2]=1;temp[3]=6;
> > >   ierr = MatSetValues(mat,2,ind1,2,ind1,temp,INSERT_VALUES);CHKERRQ(ierr);
> > >
> > >    ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
> > >    ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
> > >
> > >    MatView(mat,PETSC_VIEWER_STDOUT_SELF);
> > >    // begin cholesky factorization
> > >    IS perm,colp;
> > >
> > >    MatGetOrdering(mat,MATORDERING_NATURAL,&perm,&colp);
> > >    ierr = ISDestroy(colp);CHKERRQ(ierr);
> > >
> > >    MatFactorInfo info;
> > >    info.fill=1.0;
> > >    Mat result;
> > >
> > >    ierr = MatCholeskyFactorSymbolic(mat,perm,&info,&result); CHKERRQ(ierr);
> > >    ierr = MatCholeskyFactorNumeric(mat,&info,&result);CHKERRQ(ierr);
> > >     MatView(result, PETSC_VIEWER_STDOUT_SELF);
> > >
> > >    ierr = ISDestroy(perm);CHKERRQ(ierr);
> > >    ierr= MatDestroy(mat);CHKERRQ(ierr);
> > >    PetscFinalize();
> > >    return 0;
> > > }
> > > ***************************************************************************
> > >
> > > the result given below:
> > > mat:
> > > row0: (0,3),(1,2),(4,1),(5,1)
> > > row1: (0,2),(1,3),(4,2),(5,1)
> > > row2: (2,4),(3,1)
> > > row3: (2,1),(3,5)
> > > row4: (4,5),(5,1)
> > > row5: (4,1),(5,6)
> > > result:
> > > row0: (4,0.2),(5,-0.2)
> > > row1: (4,-0.8) (5, -0.2)
> > > row2:
> > > row3:
> > > row4:
> > > row5:
> > > ***************************************************************************
> > >
> > >    The matrix is in SBAIJ format, block size is 2, which cholesky
> > factorization
> > > function should i call? really need that for the project! Thank you!
> > >
> > > Jun
> > >
> > >
> >
> >
>
>
>
>