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 <[email protected]>: > 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/8f21f52c465b775a76cda90fe9c51d > 0c742078c7 > > 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 >>>> > >>>> >>>> >>> > -- Stefano
