<div dir="ltr">I think we have a fix for this. You can test by using the branch. In PETSC_DIR:<div><br></div><div>$ git checkout mark/ksp-zero-eig<br></div><div>$ make</div><div><br></div><div>Now remake your code and test,</div><div><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Wed, Oct 22, 2014 at 8:36 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:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Granted that I solved my problem by **not** solving the system if the Rhs 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>
for(int i = 1; i < N; ++i) {<br>
to<br>
for(int i = 0; i < N; ++i) {<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 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 -ksp_monitor -<br>
ksp_monitor_true_residual<br>
0 KSP Residual norm 0.000000000000e+00<br>
0 KSP preconditioned resid norm 0.000000000000e+00 true resid norm<br>
0.000000000000e+00 ||r(i)||/||b|| -nan<br>
0 KSP Residual norm 3.964678029172e+02<br>
0 KSP preconditioned resid norm 3.964678029172e+02 true resid norm<br>
1.118033988750e+00 ||r(i)||/||b|| 1.000000000000e+00<br>
1 KSP Residual norm 6.290716828453e-13<br>
1 KSP preconditioned resid norm 6.290716828453e-13 true resid norm<br>
3.788861141557e+00 ||r(i)||/||b|| 3.388860427931e+00<br>
0 KSP Residual norm 7.929356058344e+02<br>
0 KSP preconditioned resid norm 7.929356058344e+02 true resid norm<br>
2.236067977500e+00 ||r(i)||/||b|| 1.000000000000e+00<br>
1 KSP Residual norm 1.258143365691e-12<br>
1 KSP preconditioned resid norm 1.258143365691e-12 true resid norm<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 -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> 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 --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 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 -ksp_monitor_true_residual<br>
0 KSP Residual norm 0.000000000000e+00<br>
0 KSP preconditioned resid norm 0.000000000000e+00 true resid norm<br>
0.000000000000e+00 ||r(i)||/||b|| -nan<br>
0 KSP Residual norm 1.839181443725e+00<br>
0 KSP preconditioned resid norm 1.839181443725e+00 true resid norm<br>
1.118033988750e+00 ||r(i)||/||b|| 1.000000000000e+00<br>
1 KSP Residual norm 1.767908128952e+00<br>
1 KSP preconditioned resid norm 1.767908128952e+00 true resid norm<br>
1.182712814433e+00 ||r(i)||/||b|| 1.057850500373e+00<br>
[...]<br>
82 KSP Residual norm 1.642117972681e-05<br>
82 KSP preconditioned resid norm 1.642117972681e-05 true resid norm<br>
1.556660502826e-05 ||r(i)||/||b|| 1.392319480883e-05<br>
[...]<br>
0 KSP Residual norm 3.678362887450e+00<br>
0 KSP preconditioned resid norm 3.678362887450e+00 true resid norm<br>
2.236067977500e+00 ||r(i)||/||b|| 1.000000000000e+00<br>
1 KSP Residual norm 3.535816257904e+00<br>
1 KSP preconditioned resid norm 3.535816257904e+00 true resid norm<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 ctx.rhs,<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 I<br>
> > never override ctx.rhs.<br>
> Oh yes. I missed that.<br>
><br>
> What I observed was when I ran your code each solve had an initial<br>
> residual norm of exactly zero. I’m not sure why that happens.<br>
> > For 2): that was a small typo, problem is present no matter what number of<br>
> > levels I use.<br>
> ><br>
> > Thanks,<br>
> > Filippo<br>
<div class="HOEnZb"><div class="h5">> ><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 called<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 hand<br>
> >> side.<br>
> >><br>
> >> 2) -pc_type mg doesn’t give enough information for the multigrid to know<br>
> >> how many levels to use. In fact it is always using one hence it is not<br>
> >> multigrid. You need to use -p_mg_levels 2 or 3 or whatever to get it to<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 and<br>
> >>>>> multiple solves. I got a linear system that I solve many times, with<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 MG<br>
> >>>>> works if I destroy the ksp each time or I recompute the matrices 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); CHKERRQ(ierr);<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 gamg)<br>
> >>>>> everything works (So the matrix and RHS are working fine).<br>
> >>>>> - But the problem persists with Galerkin matrices (-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 solution<br>
> >>>>> is<br>
> >>>>> also perfect (So is not a boundary scaling or other stuff problem).<br>
> >>>>> - If I run MG with only 2 levels (so just coarse) I also get a fine<br>
> >>>>> result.<br>
> >>>>><br>
> >>>>> Setup RHS is called each time as expected, setup matrix is called just<br>
> >>>>> once, also as expected.<br>
> >>>>><br>
> >>>>> The only thing I can think is that MG does not update some value that<br>
> >>>>> actually needs to be recomputed.<br>
> >>>>><br>
> >>>>> Any idea?<br>
> >>>>><br>
> >>>>> The solution is not that different from:<br>
</div></div><div class="HOEnZb"><div class="h5">> >>>>> <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>
> >>>>> als/<br>
> >>>>> ex29.c.html<br>
> >>>>><br>
> >>>>> Best,<br>
> >>>>> Filippo<br>
> >>><br>
> >>> <test1.cpp><br>
<br>
</div></div><span class="HOEnZb"><font color="#888888">--<br>
Filippo Leonardi<br>
SAM ETH Zürich<br>
(HG J 45) Rämistrasse 101, Zürich (CH)</font></span></blockquote></div><br></div>