<html><head><meta http-equiv="Content-Type" content="text/html charset=utf-8"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space;" class=""><meta http-equiv="Content-Type" content="text/html charset=utf-8" class=""><div style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space;" class=""><br class=""><div class=""><blockquote type="cite" class=""><div class="">On Sep 29, 2017, at 6:19 PM, Hong <<a href="mailto:hzhang@mcs.anl.gov" class="">hzhang@mcs.anl.gov</a>> wrote:</div><br class="Apple-interchange-newline"><div class=""><div dir="ltr" class="">Thanks for all the input. I can do following:<div class="">1) Move test to <span style="font-size:12.8px" class="">MatGetFactor()</span><span style="font-size:12.8px" class=""> </span></div><div class=""><span style="font-size:12.8px" class=""> - If there is sufficient requests from user, we are able to add Hermitian support to petsc sequential Cholesky. </span></div><div class=""><span style="font-size:12.8px" class=""><br class=""></span></div><div class=""><span style="font-size:12.8px" class="">2) Tests</span><span style="font-size:12.8px" class="">:</span></div><div class=""><span style="font-size:12.8px" class="">a) complex build && </span><span style="font-size:12.8px" class="">ftype==</span><span style="font-size:12.8px" class="">CHOLESKY</span><span style="font-size:12.8px" class="">:</span></div><div class=""><span style="font-size:12.8px" class=""> if (mat->hermitian) </span><span style="font-size:12.8px" class="">throw an "not supported" error</span></div><div class=""><span style="font-size:12.8px" class=""> </span></div><div class=""><span style="font-size:12.8px" class="">b) all builds:</span></div><div class=""><span style="font-size:12.8px" class=""> if (!sbaij && (</span><span style="font-size:12.8px" class="">CHOLESKY || ICC</span><span style="font-size:12.8px" class="">))</span></div><div class=""><span style="font-size:12.8px" class=""> if (!mat->symmetric) </span></div></div></div></blockquote><div class=""><br class=""></div><div class="">I guess checking for mat->symmetric is redundant as you are calling MatIsSymmetric that already checks for it.</div><div class="">Also, I’d rather prefer to guard this code branch with PETSC_USE_DEBUG</div><div class=""><br class=""></div><blockquote type="cite" class=""><div class=""><div dir="ltr" class=""><div class=""><span style="font-size:12.8px" class=""> call MatIsSymmetric(mat,0.0,&symm) </span></div><div class=""><span style="font-size:12.8px" class=""> if (!symm) </span><span style="font-size:12.8px" class="">throw an error</span></div><div class=""><span style="font-size:12.8px" class=""><br class=""></span></div><div class=""><span style="font-size:12.8px" class="">Hong</span></div><div class=""><span style="font-size:12.8px" class=""><br class=""></span></div></div><div class="gmail_extra"><br class=""><div class="gmail_quote">On Fri, Sep 29, 2017 at 9:55 AM, Greg Meyer <span dir="ltr" class=""><<a href="mailto:gregory.meyer@gmail.com" target="_blank" class="">gregory.meyer@gmail.com</a>></span> wrote:<br class=""><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><p dir="ltr" class="">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. </p><div class="HOEnZb"><div class="h5">
<br class=""><div class="gmail_quote"><div dir="ltr" class="">On Fri, Sep 29, 2017, 7:37 AM Barry Smith <<a href="mailto:bsmith@mcs.anl.gov" target="_blank" class="">bsmith@mcs.anl.gov</a>> wrote:<br class=""></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><br class="">
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.<br class="">
<br class="">
2) It is a tough call if we should do this check or not. There are good arguments in both directions.<br class="">
<br class="">
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.<br class="">
<br class="">
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).<br class="">
<br class="">
Hong can you move the test up to the MatGetFactor() level?<br class="">
<br class="">
Thanks<br class="">
Barry<br class="">
<br class="">
<br class="">
> On Sep 28, 2017, at 11:35 PM, Stefano Zampini <<a href="mailto:stefano.zampini@gmail.com" target="_blank" class="">stefano.zampini@gmail.com</a>> wrote:<br class="">
><br class="">
> Hong,<br class="">
><br class="">
> I personally believe that commit <a href="https://bitbucket.org/petsc/petsc/commits/966c94c9cf4fa13d455726ec36800a577e00b171" rel="noreferrer" target="_blank" class="">https://bitbucket.org/petsc/<wbr class="">petsc/commits/<wbr class="">966c94c9cf4fa13d455726ec36800a<wbr class="">577e00b171</a> should be reverted.<br class="">
> 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 (<a href="https://bitbucket.org/petsc/petsc/commits/8f21f52c465b775a76cda90fe9c51d0c742078c7" rel="noreferrer" target="_blank" class="">https://bitbucket.org/petsc/<wbr class="">petsc/commits/<wbr class="">8f21f52c465b775a76cda90fe9c51d<wbr class="">0c742078c7</a>) , 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 <a href="https://bitbucket.org/petsc/petsc/commits/966c94c9cf4fa13d455726ec36800a577e00b171" rel="noreferrer" target="_blank" class="">https://bitbucket.org/petsc/<wbr class="">petsc/commits/<wbr class="">966c94c9cf4fa13d455726ec36800a<wbr class="">577e00b171</a>).<br class="">
><br class="">
> Barry, what do you think?<br class="">
><br class="">
> 2017-09-29 5:28 GMT+03:00 Hong <<a href="mailto:hzhang@mcs.anl.gov" target="_blank" class="">hzhang@mcs.anl.gov</a>>:<br class="">
> Greg:<br class="">
> 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! :)<br class="">
><br class="">
> I added an error flag for using MUMPS Cholesky factorization on Hermitian matrix. See branch hzhang/mumps-<wbr class="">HermitianCholesky-errflag<br class="">
> <a href="https://bitbucket.org/petsc/petsc/commits/8f21f52c465b775a76cda90fe9c51d0c742078c7" rel="noreferrer" target="_blank" class="">https://bitbucket.org/petsc/<wbr class="">petsc/commits/<wbr class="">8f21f52c465b775a76cda90fe9c51d<wbr class="">0c742078c7</a><br class="">
><br class="">
> Jose,<br class="">
> PETSc does not support Cholesky for Hermitian matrix.<br class="">
><br class="">
> 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/<wbr class="">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:<br class="">
> ierr = MatSetValue(A,0,1,1.0+PETSC_i,<wbr class="">INSERT_VALUES);CHKERRQ(ierr);<br class="">
> ierr = MatSetValue(A,1,0,1.0-PETSC_i,<wbr class="">INSERT_VALUES);CHKERRQ(ierr);<br class="">
><br class="">
> In this case, you must call<br class="">
> MatSetOption(A,MAT_HERMITIAN,<wbr class="">PETSC_TRUE);<br class="">
><br class="">
> Then, petsc will throw an error for '-pc_type cholesky'.<br class="">
><br class="">
> 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.<br class="">
><br class="">
> I also add an error flag for Cholesky/ICC if user does not set<br class="">
> MatSetOption(A,MAT_SYMMETRIC,<wbr class="">PETSC_TRUE) for aij matrix.<br class="">
> See <a href="https://bitbucket.org/petsc/petsc/commits/966c94c9cf4fa13d455726ec36800a577e00b171" rel="noreferrer" target="_blank" class="">https://bitbucket.org/petsc/<wbr class="">petsc/commits/<wbr class="">966c94c9cf4fa13d455726ec36800a<wbr class="">577e00b171</a><br class="">
><br class="">
> Let me know if you have any comments about this fix.<br class="">
><br class="">
> Hong<br class="">
><br class="">
> > El 25 sept 2017, a las 7:17, Greg Meyer <<a href="mailto:gregory.meyer@gmail.com" target="_blank" class="">gregory.meyer@gmail.com</a>> escribió:<br class="">
> ><br class="">
> > Hi all,<br class="">
> ><br class="">
> > 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: <a href="https://www.mcs.anl.gov/petsc/documentation/linearsolvertable.html" rel="noreferrer" target="_blank" class="">https://www.mcs.anl.gov/petsc/<wbr class="">documentation/<wbr class="">linearsolvertable.html</a><br class="">
> ><br class="">
> > So do you or anyone have any idea why I get incorrect eigenvalues?<br class="">
> ><br class="">
> > Thanks,<br class="">
> > Greg<br class="">
> ><br class="">
> > On Thu, Sep 21, 2017 at 5:51 PM Greg Meyer <<a href="mailto:gregory.meyer@gmail.com" target="_blank" class="">gregory.meyer@gmail.com</a>> wrote:<br class="">
> > 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...<br class="">
> ><br class="">
> > On Thu, Sep 21, 2017 at 5:24 PM Hong <<a href="mailto:hzhang@mcs.anl.gov" target="_blank" class="">hzhang@mcs.anl.gov</a>> wrote:<br class="">
> > Greg :<br class="">
> > Yes, they are Hermitian.<br class="">
> ><br class="">
> > PETSc does not support Cholesky factorization for Hermitian.<br class="">
> > It seems mumps does not support Hermitian either<br class="">
> > <a href="https://lists.mcs.anl.gov/mailman/htdig/petsc-users/2015-November/027541.html" rel="noreferrer" target="_blank" class="">https://lists.mcs.anl.gov/<wbr class="">mailman/htdig/petsc-users/<wbr class="">2015-November/027541.html</a><br class="">
> ><br class="">
> > Hong<br class="">
> ><br class="">
> ><br class="">
> > On Thu, Sep 21, 2017 at 3:43 PM Hong <<a href="mailto:hzhang@mcs.anl.gov" target="_blank" class="">hzhang@mcs.anl.gov</a>> wrote:<br class="">
> > Greg:<br class="">
> ><br class="">
> > 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!<br class="">
> > Are your matrices symmetric/Hermitian?<br class="">
> > Hong<br class="">
> ><br class="">
> > On Thu, Sep 21, 2017 at 2:05 PM Hong <<a href="mailto:hzhang@mcs.anl.gov" target="_blank" class="">hzhang@mcs.anl.gov</a>> wrote:<br class="">
> > Gregory :<br class="">
> > Use '-eps_view' for both runs to check the algorithms being used.<br class="">
> > Hong<br class="">
> ><br class="">
> > Hi all,<br class="">
> ><br class="">
> > I'm using shift-invert with EPS to solve for eigenvalues. I find that if I do only<br class="">
> ><br class="">
> > ...<br class="">
> > ierr = EPSGetST(eps,&st);CHKERRQ(<wbr class="">ierr);<br class="">
> > ierr = STSetType(st,STSINVERT);<wbr class="">CHKERRQ(ierr);<br class="">
> > ...<br class="">
> ><br class="">
> > in my code I get correct eigenvalues. But if I do<br class="">
> ><br class="">
> > ...<br class="">
> > ierr = EPSGetST(eps,&st);CHKERRQ(<wbr class="">ierr);<br class="">
> > ierr = STSetType(st,STSINVERT);<wbr class="">CHKERRQ(ierr);<br class="">
> > ierr = STGetKSP(st,&ksp);CHKERRQ(<wbr class="">ierr);<br class="">
> > ierr = KSPGetPC(ksp,&pc);CHKERRQ(<wbr class="">ierr);<br class="">
> > ierr = KSPSetType(ksp,KSPPREONLY);<wbr class="">CHKERRQ(ierr);<br class="">
> > ierr = PCSetType(pc,PCCHOLESKY);<wbr class="">CHKERRQ(ierr);<br class="">
> > ...<br class="">
> ><br class="">
> > 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.<br class="">
> ><br class="">
> > Best,<br class="">
> > Greg<br class="">
> ><br class="">
><br class="">
><br class="">
><br class="">
><br class="">
><br class="">
><br class="">
> --<br class="">
> Stefano<br class="">
<br class="">
</blockquote></div>
</div></div></blockquote></div><br class=""></div>
</div></blockquote></div><br class=""></div></body></html>