[petsc-dev] Complex Cholesky

[Jose E. Roman]

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