[petsc-dev] Complex Cholesky

[Hong Zhang]

Jose,
sbaij is for symmetric matrix, not Hermitian, as you noticed.

Petsc does not have support for Hermitian operations for Cholesky solve yet.
We need MatSolve_SeqSBAIJ_1_Hermitian().

Hong

On Wed, Jul 25, 2012 at 3:08 PM, Jose E. Roman <jroman at dsic.upv.es> wrote:

> In complex arithmetic, Cholesky requires a complex Hermitian positive
> definite matrix, right?
>
> $ ./ex10 -f
> $PETSC_DIR/share/petsc/datafiles/matrices/hpd-complex-int32-float64
> -ksp_type preonly -pc_type lu
> Number of iterations =   1
> Residual norm < 1.e-12
>
> $ ./ex10 -f
> $PETSC_DIR/share/petsc/datafiles/matrices/hpd-complex-int32-float64
> -ksp_type preonly -pc_type cholesky
> Number of iterations =   1
> Residual norm 1.64361
>
> The second run gives a wrong solution. I guess I should run with SBAIJ:
>
> $ ./ex10 -f
> $PETSC_DIR/share/petsc/datafiles/matrices/hpd-complex-int32-float64
> -ksp_type preonly -pc_type cholesky -mat_type sbaij
> Number of iterations =   1
> Residual norm < 1.e-12
>
> Ok. The thing is that this latter run is not solving the system that I
> intended to, i.e., the matrix is not the same as in the first run. If I add
> -mat_view_matlab I can see that the matrix is not complex Hermitian but
> complex symmetric, because SBAIJ only reads the upper triangular part and
> assumes that the matrix is symmetric. If I add
> MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE) after MatLoad as hinted in the
> MATSEQSBAIJ manpage, I get a bad residual:
>
> $ ./ex10 -f
> $PETSC_DIR/share/petsc/datafiles/matrices/hpd-complex-int32-float64
> -ksp_type preonly -pc_type cholesky -mat_type sbaij
> Number of iterations =   1
> Residual norm 1.64361
>
> ... and -mat_view_matlab does not show a complex Hermitian matrix either.
>
> I am a bit confused. Am I doing something wrong? Or is there a bug
> somewhere?
>
> Jose
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20120725/1d9f73a4/attachment.html>