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