Should we add a basic fallback that tests for matmult and matmulttranspose to MatIsSymmetric if ops->issymmetric is not set?
Il 29 Set 2017 7:47 PM, "Hong" <[email protected]> ha scritto: > Barry, > > When users apply Cholesky to a non-symmetric matrix, > petsc uses his upper half matrix and would produce incorrect solutions > without user's knowledge. > > Adding such check under "#PETSC_USE_DEBUG", user sees > 1) his code crashes when matrix is non-symmetric > or > 2) too slow due to checking, which prompts user to set symmetric flag. > I do not see any harm to call MatIsSymmetric() in this situation. > > Hong > > On Fri, Sep 29, 2017 at 10:57 AM, Barry Smith <[email protected]> wrote: > >> >> 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/pe >> tsc/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/p >> etsc/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/pe >> tsc/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/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- >> 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 >> > > >> > > >> > >> > >> >> >
