Cholesky factorization problem

[Hong Zhang]

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
>
>