<div dir="ltr"><br><div class="gmail_extra"><br><div class="gmail_quote">On Thu, Oct 23, 2014 at 2:32 AM, Filippo Leonardi <span dir="ltr"><<a href="mailto:filippo.leonardi@sam.math.ethz.ch" target="_blank">filippo.leonardi@sam.math.ethz.ch</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">Tested on this branch: the program no longer crashes, but the solution is<br>
actually still wrong.<br>
<br>
To reproduce you can compile a version with the for cycle starting from i = 0<br>
and one from i = 1 and see the difference with:</blockquote><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
<span class="">mpirun -np 4 ./test1 -pc_type mg -pc_mg_levels 5 -draw_pause -1<br>
</span>and without mg.<br>
<div class=""><div class="h5"><br></div></div></blockquote><div><br></div><div>I'm not sure what you are asking. starting i=0 or i=1, and with and without mg. So four possibilities? The solutions will be different for i=0 or 1 in this code:</div><div><br></div><div><div><div> ierr = VecSet(ctx.rhs, 1); CHKERRQ(ierr);</div><div> </div><div> ierr = setup(&ctx); CHKERRQ(ierr);</div><div> </div><div> PetscViewer viewer;</div><div> PetscViewerDrawOpen(PETSC_COMM_WORLD,NULL,NULL,0,0,300,300,&viewer);</div><div> for(int i = 0; i < N; ++i) {</div><div> </div><div> ierr = VecSet(ctx.rhs, (float) i * mx); CHKERRQ(ierr);</div><div> ierr = VecScale(ctx.rhs, rank); CHKERRQ(ierr);</div><div> ierr = KSPSolve(ksp, NULL, NULL); CHKERRQ(ierr);</div><div> ierr = KSPGetSolution(ksp, &phi); CHKERRQ(ierr);</div><div> ierr = VecView(ctx.rhs, viewer); CHKERRQ(ierr);</div><div> ierr = VecView(phi, viewer); CHKERRQ(ierr);</div><div> }</div></div><div><br></div><div> </div></div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div class=""><div class="h5">
On Wednesday 22 October 2014 09:31:40 Mark Adams wrote:<br>
> I think we have a fix for this. You can test by using the branch. In<br>
> PETSC_DIR:<br>
><br>
> $ git checkout mark/ksp-zero-eig<br>
> $ make<br>
><br>
> Now remake your code and test,<br>
><br>
><br>
> On Wed, Oct 22, 2014 at 8:36 AM, Filippo Leonardi <<br>
><br>
> <a href="mailto:filippo.leonardi@sam.math.ethz.ch">filippo.leonardi@sam.math.ethz.ch</a>> wrote:<br>
> > Granted that I solved my problem by **not** solving the system if the Rhs<br>
> > is<br>
> > zero. Here's some result:<br>
> ><br>
> > 1) Using Multigrid with debug crashes (see below)<br>
> > 2) Using MG without debug doesn't crash, but gives wrong solution<br>
> > 3) defaults or gamg or wathever*works both with and without debug<br>
> ><br>
> > Now weirder:<br>
> > 4) changing line 145 from<br>
> ><br>
> > for(int i = 1; i < N; ++i) {<br>
> ><br>
> > to<br>
> ><br>
> > for(int i = 0; i < N; ++i) {<br>
> ><br>
> > i.e. avoiding broken case 1) above (rhs never 0) makes 1), 2) and 3) work<br>
> > perfectly (and any combination thereof)<br>
> ><br>
> > (You can look at the differences in the pictures given by Draw to see that<br>
> > the<br>
> > solution is more or less wrong.)<br>
> ><br>
> > Details:<br>
> ><br>
> > Release build with no debug:<br>
> ><br>
> > mpirun -np 4 ./test1 -pc_type mg -pc_mg_levels 5 -draw_pause -1<br>
> > -ksp_monitor -<br>
> > ksp_monitor_true_residual<br>
> ><br>
> > 0 KSP Residual norm 0.000000000000e+00<br>
> > 0 KSP preconditioned resid norm 0.000000000000e+00 true resid norm<br>
> ><br>
> > 0.000000000000e+00 ||r(i)||/||b|| -nan<br>
> ><br>
> > 0 KSP Residual norm 3.964678029172e+02<br>
> > 0 KSP preconditioned resid norm 3.964678029172e+02 true resid norm<br>
> ><br>
> > 1.118033988750e+00 ||r(i)||/||b|| 1.000000000000e+00<br>
> ><br>
> > 1 KSP Residual norm 6.290716828453e-13<br>
> > 1 KSP preconditioned resid norm 6.290716828453e-13 true resid norm<br>
> ><br>
> > 3.788861141557e+00 ||r(i)||/||b|| 3.388860427931e+00<br>
> ><br>
> > 0 KSP Residual norm 7.929356058344e+02<br>
> > 0 KSP preconditioned resid norm 7.929356058344e+02 true resid norm<br>
> ><br>
> > 2.236067977500e+00 ||r(i)||/||b|| 1.000000000000e+00<br>
> ><br>
> > 1 KSP Residual norm 1.258143365691e-12<br>
> > 1 KSP preconditioned resid norm 1.258143365691e-12 true resid norm<br>
> ><br>
> > 7.577722283114e+00 ||r(i)||/||b|| 3.388860427931e+00<br>
> ><br>
> > i.e. I don't see the residual to be zero each time.<br>
> ><br>
> > Debug build:<br>
> ><br>
> > mpirun -np 4 ./test1 -pc_type mg -pc_mg_levels 5 -draw_pause -1<br>
> > -ksp_monitor -<br>
> > ksp_monitor_true_residual<br>
> > [1]PETSC ERROR: [0]PETSC ERROR: --------------------- Error Message<br>
> > --------------------------------------------------------------<br>
> > [0]PETSC ERROR: Null argument, when expecting valid pointer<br>
> > [0]PETSC ERROR: Null Pointer: Parameter # 3<br>
> > [0]PETSC ERROR: See <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html" target="_blank">http://www.mcs.anl.gov/petsc/documentation/faq.html</a><br>
> > for<br>
> > trouble shooting.<br>
> > [0]PETSC ERROR: Petsc Release Version 3.5.0, Jun, 30, 2014<br>
> > [0]PETSC ERROR: ./test1 on a arch-linux2-cxx-debug named Besikovitch-II by<br>
> > filippo Wed Oct 22 14:20:46 2014<br>
> > [0]PETSC ERROR: Configure options --with-clanguage=cxx<br>
> > --download-scalapack --<br>
> > download-mumps --download-parmetis --download-metis --download-fftw<br>
> > [0]PETSC ERROR: #1 KSPComputeEigenvalues() line 119 in<br>
> > /home/filippo/Workspace/petsc-3.5.0/src/ksp/ksp/interface/itfunc.c<br>
> > [0]PETSC ERROR: #2 KSPChebyshevComputeExtremeEigenvalues_Private() line<br>
> > 328 in<br>
> > /home/filippo/Workspace/petsc-3.5.0/src/ksp/ksp/impls/cheby/cheby.c<br>
> > [0]PETSC ERROR: #3 KSPSolve_Chebyshev() line 373 in<br>
> > /home/filippo/Workspace/petsc-3.5.0/src/ksp/ksp/impls/cheby/cheby.c<br>
> > [0]PETSC ERROR: #4 KSPSolve() line 459 in<br>
> > /home/filippo/Workspace/petsc-3.5.0/src/ksp/ksp/interface/itfunc.c<br>
> > [0]PETSC ERROR: #5 PCMGMCycle_Private() line 19 in<br>
> > /home/filippo/Workspace/petsc-3.5.0/src/ksp/pc/impls/mg/mg.c<br>
> ><br>
> > Debug build again with defaults:<br>
> ><br>
> > mpirun -np 4 ./test1 -draw_pause -1 -ksp_monitor<br>
> > -ksp_monitor_true_residual<br>
> ><br>
> > 0 KSP Residual norm 0.000000000000e+00<br>
> > 0 KSP preconditioned resid norm 0.000000000000e+00 true resid norm<br>
> ><br>
> > 0.000000000000e+00 ||r(i)||/||b|| -nan<br>
> ><br>
> > 0 KSP Residual norm 1.839181443725e+00<br>
> > 0 KSP preconditioned resid norm 1.839181443725e+00 true resid norm<br>
> ><br>
> > 1.118033988750e+00 ||r(i)||/||b|| 1.000000000000e+00<br>
> ><br>
> > 1 KSP Residual norm 1.767908128952e+00<br>
> > 1 KSP preconditioned resid norm 1.767908128952e+00 true resid norm<br>
> ><br>
> > 1.182712814433e+00 ||r(i)||/||b|| 1.057850500373e+00<br>
> > [...]<br>
> ><br>
> > 82 KSP Residual norm 1.642117972681e-05<br>
> > 82 KSP preconditioned resid norm 1.642117972681e-05 true resid norm<br>
> ><br>
> > 1.556660502826e-05 ||r(i)||/||b|| 1.392319480883e-05<br>
> > [...]<br>
> > 0 KSP Residual norm 3.678362887450e+00<br>
> ><br>
> > 0 KSP preconditioned resid norm 3.678362887450e+00 true resid norm<br>
> ><br>
> > 2.236067977500e+00 ||r(i)||/||b|| 1.000000000000e+00<br>
> ><br>
> > 1 KSP Residual norm 3.535816257904e+00<br>
> > 1 KSP preconditioned resid norm 3.535816257904e+00 true resid norm<br>
> ><br>
> > 2.365425628867e+00 ||r(i)||/||b|| 1.057850500373e+00<br>
> ><br>
> > On Wednesday 22 October 2014 07:12:50 you wrote:<br>
> > > > On Oct 22, 2014, at 1:09 AM, Leonardi Filippo<br>
> > > > <<a href="mailto:filippo.leonardi@sam.math.ethz.ch">filippo.leonardi@sam.math.ethz.ch</a>> wrote:<br>
> > > ><br>
> > > > For 1): I must be missing something then: in main i set and scale<br>
> ><br>
> > ctx.rhs,<br>
> ><br>
> > > > then I call KSPSolve that calls ComputeRhs that further scales ctx.rhs<br>
> > > > and copy it to b, which is then used by the solve. To my understanding<br>
> ><br>
> > I<br>
> ><br>
> > > > never override ctx.rhs.<br>
> > > ><br>
> > > Oh yes. I missed that.<br>
> > ><br>
> > > What I observed was when I ran your code each solve had an initial<br>
> > ><br>
> > > residual norm of exactly zero. I’m not sure why that happens.<br>
> > ><br>
> > > > For 2): that was a small typo, problem is present no matter what<br>
> ><br>
> > number of<br>
> ><br>
> > > > levels I use.<br>
> > > ><br>
> > > > Thanks,<br>
> > > > Filippo<br>
> > > ><br>
> > > >> On 22-ott-2014, at 04:27, "Barry Smith" <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>> wrote:<br>
> > > >><br>
> > > >><br>
> > > >> There are a few issues with your code<br>
> > > >><br>
> > > >> 1) If you use KSPSetComputeRHS() then the function you provide is<br>
> ><br>
> > called<br>
> ><br>
> > > >> immediately before the solve automatically for every solve. Thus your<br>
> > > >> code with the vector operations here<br>
> > > >><br>
> > > >> ierr = VecSet(ctx.rhs, (float) i * mx); CHKERRQ(ierr);<br>
> > > >><br>
> > > >> ierr = VecScale(ctx.rhs, rank); CHKERRQ(ierr);<br>
> > > >> ierr = KSPSolve(ksp, NULL, NULL); CHKERRQ(ierr);<br>
> > > >><br>
> > > >> does nothing. And each linear solve ends up being with a zero right<br>
> ><br>
> > hand<br>
> ><br>
> > > >> side.<br>
> > > >><br>
> > > >> 2) -pc_type mg doesn’t give enough information for the multigrid to<br>
> ><br>
> > know<br>
> ><br>
> > > >> how many levels to use. In fact it is always using one hence it is<br>
> > > >> not<br>
> > > >> multigrid. You need to use -p_mg_levels 2 or 3 or whatever to get it<br>
> ><br>
> > to<br>
> ><br>
> > > >> use more than one level.<br>
> > > >><br>
> > > >> Barry<br>
> > > >><br>
> > > >>> On Oct 21, 2014, at 1:50 AM, Filippo Leonardi<br>
> > > >>> <<a href="mailto:filippo.leonardi@sam.math.ethz.ch">filippo.leonardi@sam.math.ethz.ch</a>> wrote:<br>
> > > >>><br>
> > > >>> Here,<br>
> > > >>><br>
> > > >>> run this e.g. with<br>
> > > >>> mpirun -np 4 ./test1 -draw_pause -1 -pc_type mg -pc_mg_galerkin<br>
> > > >>><br>
> > > >>> In Debug mode it crashes, in standard mode it fails. Using Petsc<br>
> > > >>> current.<br>
> > > >>><br>
> > > >>> Run without mg and works!<br>
> > > >>><br>
> > > >>> Hope it's not me being stupid and doing something terribly wrong!<br>
> > > >>><br>
> > > >>> Best,<br>
> > > >>> Filippo<br>
> > > >>><br>
> > > >>>> On Monday 20 October 2014 07:29:04 Barry Smith wrote:<br>
> > > >>>> Can you send us the code: to <a href="mailto:petsc-maint@mcs.anl.gov">petsc-maint@mcs.anl.gov</a> or<br>
> > > >>>> <a href="mailto:petsc-users@mcs.anl.gov">petsc-users@mcs.anl.gov</a> ? Or something that reproduces the problem?<br>
> > > >>>><br>
> > > >>>> Barry<br>
> > > >>>><br>
> > > >>>>> On Oct 20, 2014, at 3:31 AM, Filippo Leonardi<br>
> > > >>>>> <<a href="mailto:filippo.leonardi@sam.math.ethz.ch">filippo.leonardi@sam.math.ethz.ch</a>> wrote:<br>
> > > >>>>><br>
> > > >>>>> Hi,<br>
> > > >>>>><br>
> > > >>>>> I have a very specific problem that I cannot figure out with PCMG<br>
> ><br>
> > and<br>
> ><br>
> > > >>>>> multiple solves. I got a linear system that I solve many times,<br>
> ><br>
> > with<br>
> ><br>
> > > >>>>> same<br>
> > > >>>>> matrix but different RHS. I can successfully solve the system with<br>
> > > >>>>> standard techniques, such as default solver or LU or PCGAMG. Even<br>
> ><br>
> > MG<br>
> ><br>
> > > >>>>> works if I destroy the ksp each time or I recompute the matrices<br>
> > > >>>>> at<br>
> > > >>>>> each<br>
> > > >>>>> time. But when I try and go to MG and not recomputing the matrices<br>
> > > >>>>> each<br>
> > > >>>>> time the solver fails. Any idea?<br>
> > > >>>>><br>
> > > >>>>> Here some detail: the setup:<br>
> > > >>>>> ierr = KSPSetDMActive(ksp, PETSC_FALSE); CHKERRQ(ierr);<br>
> > > >>>>> ierr = KSPSetDM(ksp, da); CHKERRQ(ierr);<br>
> > > >>>>> ierr = KSPSetComputeOperators(ksp, ComputeMatrix, ctx);<br>
> ><br>
> > CHKERRQ(ierr);<br>
> ><br>
> > > >>>>> ierr = KSPSetComputeRHS(ksp, ComputeRHS, ctx); CHKERRQ(ierr);<br>
> > > >>>>> ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr);<br>
> > > >>>>><br>
> > > >>>>> Then I solve as usual, for a large number of time steps:<br>
> > > >>>>> ierr = KSPSolve(ksp, NULL, NULL); CHKERRQ(ierr);<br>
> > > >>>>> ierr = KSPGetSolution(ksp, &phi);<br>
> > > >>>>><br>
> > > >>>>> The solver converges and does that in a reasonable number of<br>
> > > >>>>> iterations:<br>
> > > >>>>> Linear solve converged due to CONVERGED_RTOL iterations 7<br>
> > > >>>>> And ksp_view and ksp_monitor do not show any weird stuff.<br>
> > > >>>>><br>
> > > >>>>> - Weirdly enough using any solver (for instance cg+bjacobi or<br>
> > > >>>>> gamg)<br>
> > > >>>>> everything works (So the matrix and RHS are working fine).<br>
> > > >>>>> - But the problem persists with Galerkin matrices<br>
> > > >>>>> (-pc_mg_galerkin)<br>
> > > >>>>> (So is<br>
> > > >>>>> not a ComputeMatrix problem).<br>
> > > >>>>> - If I do:<br>
> > > >>>>> ierr = KSPSetComputeOperators(ksp, ComputeMatrix, this);<br>
> > > >>>>> CHKERRQ(ierr);<br>
> > > >>>>> between each solve or destroy the ksp entirely each time the<br>
> ><br>
> > solution<br>
> ><br>
> > > >>>>> is<br>
> > > >>>>> also perfect (So is not a boundary scaling or other stuff<br>
> > > >>>>> problem).<br>
> > > >>>>> - If I run MG with only 2 levels (so just coarse) I also get a<br>
> > > >>>>> fine<br>
> > > >>>>> result.<br>
> > > >>>>><br>
> > > >>>>> Setup RHS is called each time as expected, setup matrix is called<br>
> ><br>
> > just<br>
> ><br>
> > > >>>>> once, also as expected.<br>
> > > >>>>><br>
> > > >>>>> The only thing I can think is that MG does not update some value<br>
> ><br>
> > that<br>
> ><br>
> > > >>>>> actually needs to be recomputed.<br>
> > > >>>>><br>
> > > >>>>> Any idea?<br>
> ><br>
> > > >>>>> The solution is not that different from:<br>
> > <a href="http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutori" target="_blank">http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutori</a><br>
> ><br>
> > > >>>>> als/<br>
> > > >>>>> ex29.c.html<br>
> > > >>>>><br>
> > > >>>>> Best,<br>
> > > >>>>> Filippo<br>
> > > >>><br>
> > > >>> <test1.cpp><br>
> ><br>
> > --<br>
> > Filippo Leonardi<br>
> > SAM ETH Zürich<br>
> > (HG J 45) Rämistrasse 101, Zürich (CH)<br>
<br>
--<br>
Filippo Leonardi<br>
SAM ETH Zürich<br>
(HG J 45) Rämistrasse 101, Zürich (CH)</div></div></blockquote></div><br></div></div>