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

[Stefano Zampini]

Perhaps you can check for the hermitian option and complex values
in MatGetFactor_XXX_mumps?

https://bitbucket.org/petsc/petsc/src/8f21f52c465b775a76cda90fe9c51d0c742078c7/src/mat/impls/aij/mpi/mumps/mumps.c?at=hzhang%2Fmumps-HermitianCholesky-errflag&fileviewer=file-view-default#mumps.c-2159

2017-09-29 9:35 GMT+03:00 Stefano Zampini <stefano.zampini at gmail.com>:

> 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/8f21f52c465b775a76cda90fe9c51d
> 0c742078c7) , 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/966c94c9cf
> 4fa13d455726ec36800a577e00b171).
>
> 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-Nov
>>>>> ember/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
>



-- 
Stefano
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170929/a4c6a5b7/attachment-0001.html>