[petsc-users] Incorrect Eigenvalues when Setting KSP and PC types

[Hong]

Greg:

> Thanks so much for the detailed response. I am glad to know the reason
> behind it--hopefully we eventually figure out why the solvers have this
> behavior! Hong, I really appreciate you working on a patch to throw an
> error in this case. It definitely bit me and some people using my code...
> Hopefully it doesn't happen to anyone else! :)
>
I added an error flag for using MUMPS Cholesky factorization on Hermitian
matrix. See branch hzhang/mumps-HermitianCholesky-errflag
https://bitbucket.org/petsc/petsc/commits/8f21f52c465b775a76cda90fe9c51d0c742078c7


Jose,
PETSc does not support Cholesky for Hermitian matrix.

>
>>> The linear solver table probably needs to be updated. I have tried
>>> several Cholesky solvers. With mkl_pardiso I get an explicit error message
>>> that it does not support Cholesky with complex scalars. The rest (PETSc,
>>> MUMPS, CHOLMOD) give a wrong answer (without error message). The problem is
>>> not related to your matrix, nor to shift-and-invert in SLEPc. I tried with
>>> ex1.c under PETSC_DIR/src/ksp/ksp/examples/tutorials. The example works
>>> in complex scalars, but the matrix is real. As soon as you add complex
>>> entries Cholesky fails, for instance adding this:
>>>   ierr = MatSetValue(A,0,1,1.0+PETSC_i,INSERT_VALUES);CHKERRQ(ierr);
>>>   ierr = MatSetValue(A,1,0,1.0-PETSC_i,INSERT_VALUES);CHKERRQ(ierr);
>>>
>>
In this case, you must call
MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE);

Then, petsc will throw an error for '-pc_type cholesky'.

>
>>> I don't know if it is a bug or that the algorithm cannot support complex
>>> Hermitian matrices. Maybe Hong can confirm any of these. In the latter
>>> case, I agree that all packages should give an error message, as
>>> mkl_pardiso does.
>>>
>>
I also add an error flag for Cholesky/ICC if user does not set
MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE) for aij matrix.
See
https://bitbucket.org/petsc/petsc/commits/966c94c9cf4fa13d455726ec36800a577e00b171

Let me know if you have any comments about this fix.

Hong

>
>>> > El 25 sept 2017, a las 7:17, Greg Meyer <gregory.meyer at gmail.com>
>>> escribió:
>>> >
>>> > Hi all,
>>> >
>>> > Hong--looking at your link, there may be no special algorithm for
>>> Hermitian matrices in MUMPS, but that doesn't mean it can't solve them like
>>> it would any matrix. Furthermore it appears that Cholesky of complex
>>> matrices is supported from this link: https://www.mcs.anl.gov/petsc/
>>> documentation/linearsolvertable.html
>>> >
>>> > So do you or anyone have any idea why I get incorrect eigenvalues?
>>> >
>>> > Thanks,
>>> > Greg
>>> >
>>> > On Thu, Sep 21, 2017 at 5:51 PM Greg Meyer <gregory.meyer at gmail.com>
>>> wrote:
>>> > Ok, thanks. It seems that PETSc clearly should throw an error in this
>>> case instead of just giving incorrect answers? I am surprised that it does
>>> not throw an error...
>>> >
>>> > On Thu, Sep 21, 2017 at 5:24 PM Hong <hzhang at mcs.anl.gov> wrote:
>>> > Greg :
>>> > Yes, they are Hermitian.
>>> >
>>> > PETSc does not support  Cholesky factorization for Hermitian.
>>> > It seems mumps does not support Hermitian either
>>> > https://lists.mcs.anl.gov/mailman/htdig/petsc-users/
>>> 2015-November/027541.html
>>> >
>>> > Hong
>>> >
>>> >
>>> > On Thu, Sep 21, 2017 at 3:43 PM Hong <hzhang at mcs.anl.gov> wrote:
>>> > Greg:
>>> >
>>> > OK, the difference is whether LU or Cholesky factorization is used.
>>> But I would hope that neither one should give incorrect eigenvalues, and
>>> when I run with the latter it does!
>>> > Are your matrices symmetric/Hermitian?
>>> > Hong
>>> >
>>> > On Thu, Sep 21, 2017 at 2:05 PM Hong <hzhang at mcs.anl.gov> wrote:
>>> > Gregory :
>>> > Use '-eps_view' for both runs to check the algorithms being used.
>>> > Hong
>>> >
>>> > Hi all,
>>> >
>>> > I'm using shift-invert with EPS to solve for eigenvalues. I find that
>>> if I do only
>>> >
>>> > ...
>>> >   ierr = EPSGetST(eps,&st);CHKERRQ(ierr);
>>> >   ierr = STSetType(st,STSINVERT);CHKERRQ(ierr);
>>> > ...
>>> >
>>> > in my code I get correct eigenvalues. But if I do
>>> >
>>> > ...
>>> >   ierr = EPSGetST(eps,&st);CHKERRQ(ierr);
>>> >   ierr = STSetType(st,STSINVERT);CHKERRQ(ierr);
>>> >   ierr = STGetKSP(st,&ksp);CHKERRQ(ierr);
>>> >   ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
>>> >   ierr = KSPSetType(ksp,KSPPREONLY);CHKERRQ(ierr);
>>> >   ierr = PCSetType(pc,PCCHOLESKY);CHKERRQ(ierr);
>>> > ...
>>> >
>>> > the eigenvalues found by EPS are completely wrong! Somehow I thought I
>>> was supposed to do the latter, from the examples etc, but I guess that was
>>> not correct? I attach the full piece of test code and a test matrix.
>>> >
>>> > Best,
>>> > Greg
>>> >
>>>
>>>
>>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170928/8e36de2c/attachment.html>