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 <[email protected]>: > 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 <[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/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 <[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-Nov >>>>> ember/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 > -- Stefano
