No, never do the MatIsSymmetric() that was me just "thinking aloud" in my first mail
Barry > On Sep 29, 2017, at 8:33 AM, Hong <[email protected]> wrote: > > Taking both suggestions, how about > > 2) Tests: > a) complex build && ftype==CHOLESKY: > if (mat->hermitian) throw an "not supported" error > > b) all builds: > if (!mat->sbaij && (CHOLESKY || ICC)) > if (!mat->symmetric) //user does not indicate mat is symmetric > #PETSC_USE_DEBUG > call MatIsSymmetric(mat,0.0,&symm) > if (!symm) throw an error > #else > throw an error > #endif > > Hong > > On Fri, Sep 29, 2017 at 10:25 AM, Barry Smith <[email protected]> wrote: > > No I don't want you to actually check if the matrix is symmetric (too > expensive) just throw an error if the user has not indicated the appropriate > properties of the matrix > > > > On Sep 29, 2017, at 8:19 AM, Hong <[email protected]> wrote: > > > > Thanks for all the input. I can do following: > > 1) Move test to MatGetFactor() > > - If there is sufficient requests from user, we are able to add > > Hermitian support to petsc sequential Cholesky. > > > > 2) Tests: > > a) complex build && ftype==CHOLESKY: > > if (mat->hermitian) throw an "not supported" error > > > > b) all builds: > > if (!sbaij && (CHOLESKY || ICC)) > > if (!mat->symmetric) > > call MatIsSymmetric(mat,0.0,&symm) > > if (!symm) throw an error > > > > Hong > > > > > > On Fri, Sep 29, 2017 at 9:55 AM, Greg Meyer <[email protected]> wrote: > > FYI my perspective as a user--something that really tricked me was that > > after setting the solver to Hermitian problem, the algorithm returned real > > eigenvalues so they seemed OK. When I turned off Hermitian as I was trying > > to debug, the eigenvalues were not at all just real, and thus it was clear > > that they were wrong. So I think the check at least when Hermitian is set > > is very important, since by construction real eigenvalues are returned. > > > > > > On Fri, Sep 29, 2017, 7:37 AM Barry Smith <[email protected]> wrote: > > > > 1) The test is definitely in the wrong place. If we are testing and > > erroring if using Cholesky and mat is not marked as symmetric or hermitian > > the test should be in MatGetFactor() not in a particular implementation. > > > > 2) It is a tough call if we should do this check or not. There are good > > arguments in both directions. > > > > One thing we could do is if the matrix is not marked as > > symmetric/hermitian is we could just check at that point (but expensive) or > > we could just check in debug mode. > > > > But I think we should require the user to set the flag (note for SBAIJ the > > flag for symmetric (or hermitian? which one) should be automatically set at > > creation). > > > > Hong can you move the test up to the MatGetFactor() level? > > > > Thanks > > Barry > > > > > > > On Sep 28, 2017, at 11:35 PM, Stefano Zampini <[email protected]> > > > wrote: > > > > > > 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/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 > > > > > > > > > > > > > > > > > > > > > > > > > -- > > > Stefano > > > > > >
