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

[Greg Meyer]

FYI my perspective as a user--something that really tricked me was that
after setting the solver to Hermitian problem, the algorithm returned real
eigenvalues so they seemed OK. When I turned off Hermitian as I was trying
to debug, the eigenvalues were not at all just real, and thus it was clear
that they were wrong. So I think the check at least when Hermitian is set
is very important, since by construction real eigenvalues are returned.

On Fri, Sep 29, 2017, 7:37 AM Barry Smith <bsmith at mcs.anl.gov> wrote:

>
>   1) The test is definitely in the wrong place. If we are testing and
> erroring if using Cholesky and mat is not marked as symmetric or hermitian
> the test should be in MatGetFactor() not in a particular implementation.
>
>  2) It is a tough call if we should do this check or not.  There are good
> arguments in both directions.
>
>    One thing we could do is if the matrix is not marked as
> symmetric/hermitian is we could just check at that point (but expensive) or
> we could just check in debug mode.
>
>  But I think we should require the user to set the flag (note for SBAIJ
> the flag for symmetric (or hermitian? which one) should be automatically
> set at creation).
>
>   Hong can you move the test up to the MatGetFactor() level?
>
>   Thanks
>     Barry
>
>
> > On Sep 28, 2017, at 11:35 PM, Stefano Zampini <stefano.zampini at gmail.com>
> wrote:
> >
> > Hong,
> >
> > I personally believe that commit
> https://bitbucket.org/petsc/petsc/commits/966c94c9cf4fa13d455726ec36800a577e00b171
> should be reverted.
> > I agree on the fact that when the user sets an option (the hermitian one
> in this case) and that feature is not supported we should throw an error (
> https://bitbucket.org/petsc/petsc/commits/8f21f52c465b775a76cda90fe9c51d0c742078c7)
> , but I don't like the fact that the user is forced to set on option to use
> a feature that he already requested (as in
> https://bitbucket.org/petsc/petsc/commits/966c94c9cf4fa13d455726ec36800a577e00b171
> ).
> >
> > Barry, what do you think?
> >
> > 2017-09-29 5:28 GMT+03:00 Hong <hzhang at mcs.anl.gov>:
> > 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
> > >
> >
> >
> >
> >
> >
> >
> > --
> > Stefano
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170929/9c2de734/attachment-0001.html>