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

Fri Sep 29 12:57:17 CDT 2017

Should we add a basic fallback that tests for matmult and matmulttranspose
to MatIsSymmetric if ops->issymmetric is not set?

Il 29 Set 2017 7:47 PM, "Hong" <hzhang at mcs.anl.gov> ha scritto:

> Barry,
>
> When users apply Cholesky to a non-symmetric matrix,
> petsc uses his upper half matrix and would produce incorrect solutions
> without user's knowledge.
>
> Adding such check under "#PETSC_USE_DEBUG", user sees
> 1) his code crashes when matrix is non-symmetric
> or
> 2) too slow due to checking, which prompts user to set symmetric flag.
> I do not see any harm to call MatIsSymmetric() in this situation.
>
> Hong
>
> On Fri, Sep 29, 2017 at 10:57 AM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
>>
>>   No, never do the MatIsSymmetric() that was me just "thinking aloud" in
>> my first mail
>>
>>   Barry
>>
>> > On Sep 29, 2017, at 8:33 AM, Hong <hzhang at mcs.anl.gov> wrote:
>> >
>> > Taking both suggestions, how about
>> >
>> > 2) Tests:
>> > a) complex build && ftype==CHOLESKY:
>> >   if (mat->hermitian) throw an "not supported" error
>> >
>> > b) all builds:
>> >   if (!mat->sbaij && (CHOLESKY || ICC))
>> >     if (!mat->symmetric) //user does not indicate mat is symmetric
>> > #PETSC_USE_DEBUG
>> >       call  MatIsSymmetric(mat,0.0,&symm)
>> >       if (!symm) throw an error
>> > #else
>> >      throw an error
>> > #endif
>> >
>> > Hong
>> >
>> > On Fri, Sep 29, 2017 at 10:25 AM, Barry Smith <bsmith at mcs.anl.gov>
>> wrote:
>> >
>> >   No I don't want you to actually check if the matrix is symmetric (too
>> expensive) just throw an error if the user has not indicated the
>> appropriate properties of the matrix
>> >
>> >
>> > > On Sep 29, 2017, at 8:19 AM, Hong <hzhang at mcs.anl.gov> wrote:
>> > >
>> > > Thanks for all the input. I can do following:
>> > > 1) Move test to MatGetFactor()
>> > >     - If there is sufficient requests from user, we are able to add
>> Hermitian support to petsc sequential Cholesky.
>> > >
>> > > 2) Tests:
>> > > a) complex build && ftype==CHOLESKY:
>> > >   if (mat->hermitian) throw an "not supported" error
>> > >
>> > > b) all builds:
>> > >   if (!sbaij && (CHOLESKY || ICC))
>> > >     if (!mat->symmetric)
>> > >       call  MatIsSymmetric(mat,0.0,&symm)
>> > >       if (!symm) throw an error
>> > >
>> > > Hong
>> > >
>> > >
>> > > On Fri, Sep 29, 2017 at 9:55 AM, Greg Meyer <gregory.meyer at gmail.com>
>> wrote:
>> > > 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/pe
>> tsc/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/p
>> etsc/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/pe
>> tsc/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/8f21f52c465b775a76
>> cda90fe9c51d0c742078c7
>> > > >
>> > > > 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/966c94c9cf4fa13d45
>> 5726ec36800a577e00b171
>> > > >
>> > > > 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/c25e3a96/attachment.html>