Hi Jose and Hong, 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! :)
Best, Greg On Mon, Sep 25, 2017, 7:44 AM Hong <[email protected]> wrote: > Greg and Jose, > I'll check this case, at least have petsc provide error message. > Please give me some time because I'm in the middle of several tasks. > I'll let you know once I add this support. > > Hong > > > On Mon, Sep 25, 2017 at 3:46 AM, Jose E. Roman <[email protected]> wrote: > >> Greg, >> >> 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); >> >> 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. >> >> As a side comment, I would suggest using LU instead of Cholesky. Cholesky >> performs less flops but it does not mean it will be faster - I have seen >> many cases where it is slower than LU, maybe because in shift-and-invert >> computations the matrix is often indefinite, so an indefinite factorization >> is computed rather than Cholesky. Some SLEPc eigensolvers (e.g. LOBPCG) >> require that the preconditioner is symmetric (Hermitian), but the default >> solver (Krylov-Schur) is quite robust when you use LU instead of Cholesky >> in Hermitian problems. And you can always solve the problem as >> non-Hermitian (the difference in accuracy should not be too noticeable). >> >> Jose >> >> >> > 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 >> > >> >> >
