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 <[email protected]> >>> 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 <[email protected]> >>> 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 <[email protected]> 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 <[email protected]> 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 <[email protected]> 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 >>> > >>> >>> >>
