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

[Greg Meyer]

Hi Jose and Hong,

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! :)

Best,
Greg

On Mon, Sep 25, 2017, 7:44 AM Hong <hzhang at mcs.anl.gov> wrote:

> Greg and Jose,
> I'll check this case, at least have petsc provide error message.
> Please give me some time because I'm in the middle of several tasks.
> I'll let you know once I add this support.
>
> Hong
>
>
> On Mon, Sep 25, 2017 at 3:46 AM, Jose E. Roman <jroman at dsic.upv.es> wrote:
>
>> Greg,
>>
>> 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);
>>
>> 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.
>>
>> As a side comment, I would suggest using LU instead of Cholesky. Cholesky
>> performs less flops but it does not mean it will be faster - I have seen
>> many cases where it is slower than LU, maybe because in shift-and-invert
>> computations the matrix is often indefinite, so an indefinite factorization
>> is computed rather than Cholesky. Some SLEPc eigensolvers (e.g. LOBPCG)
>> require that the preconditioner is symmetric (Hermitian), but the default
>> solver (Krylov-Schur) is quite robust when you use LU instead of Cholesky
>> in Hermitian problems. And you can always solve the problem as
>> non-Hermitian (the difference in accuracy should not be too noticeable).
>>
>> Jose
>>
>>
>> > 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/20170925/b33118e8/attachment.html>