[petsc-users] Matrix free GMRES seems to ignore my initial guess

Wed May 29 21:49:39 CDT 2019

Just some feedback. I found the problem. For reference my solve was called
as follows

KSPSolve(ksp,b,phi_new)

Inside my matrix operation (the "Matrix-Action" or MAT_OP_MULT) I was using
phi_new for a computation and that overwrite my initial guess everytime.
Looks like the solver still holds on to phi_new and uses it internally and
therefore when I modify it it basically changes the entire behavior. Lesson
learned: Internal to the custom MAT_OP_MULT, do not modify, b or phi_new.

Thanks for the help.

Regards,
Jan Vermaak

On Tue, May 28, 2019 at 1:16 AM Smith, Barry F. <bsmith at mcs.anl.gov> wrote:

>
>
> > On May 28, 2019, at 12:41 AM, Jan Izak Cornelius Vermaak <
> janicvermaak at gmail.com> wrote:
> >
> > Checking the matrix would be hard as I have a really big operator. Its
> transport sweeps.
> >
> > When I increase the restart interval the solution converges to the right
> one.
>
>   Run with -ksp_monitor_true_residual what are the true residuals being
> printed?
>
>   The GMRES code has been in continuous use for 25 years, it would
> stunning if you suddenly found a bug in it.
>
>    How it works, within a restart, the GMRES algorithm uses a simple
> recursive formula to compute an "estimate" for the residual norm. At
> restart it actually computes the current solution and then uses that to
> compute an accurate residual norm via the formula b - A x. When the
> residual computed by the b - A x is different than that computed by the
> recursive formula it means the recursive formula has run into some
> difficulty (bad operator, bad preconditioner, null space in the operator)
> and is not computing correct values. Now if you increase the restart to
> past the point when it "converges" you are hiding the incorrectly computed
> values computed via the recursive formula.
>
>   I urge you to check the residual norm by b - A x at the end of the solve
> and double check that it is small. It seems unlikely GMRES is providing the
> correct answer for your problem.
>
>    Barry
>
>
> > Checked against a reference solution and Classic Richardson. It is
> really as if the initial guess is completely ignored.
> >
> > [0]  Computing b
> > [0]  Iteration 0 Residual 169.302
> > [0]  Iteration 1 Residual 47.582
> > [0]  Iteration 2 Residual 13.2614
> > [0]  Iteration 3 Residual 4.46795
> > [0]  Iteration 4 Residual 1.03038
> > [0]  Iteration 5 Residual 0.246807
> > [0]  Iteration 6 Residual 0.0828341
> > [0]  Iteration 7 Residual 0.0410627
> > [0]  Iteration 8 Residual 0.0243749
> > [0]  Iteration 9 Residual 0.0136067
> > [0]  Iteration 10 Residual 0.00769078
> > [0]  Iteration 11 Residual 0.00441658
> > [0]  Iteration 12 Residual 0.00240794
> > [0]  Iteration 13 Residual 0.00132048
> > [0]  Iteration 14 Residual 0.00073003
> > [0]  Iteration 15 Residual 0.000399504
> > [0]  Iteration 16 Residual 0.000217677
> > [0]  Iteration 17 Residual 0.000120408
> > [0]  Iteration 18 Residual 6.49719e-05
> > [0]  Iteration 19 Residual 3.44523e-05
> > [0]  Iteration 20 Residual 1.87909e-05
> > [0]  Iteration 21 Residual 1.02385e-05
> > [0]  Iteration 22 Residual 5.57859e-06
> > [0]  Iteration 23 Residual 3.03431e-06
> > [0]  Iteration 24 Residual 1.63696e-06
> > [0]  Iteration 25 Residual 8.78202e-07
> >
> > On Mon, May 27, 2019 at 11:55 PM Smith, Barry F. <bsmith at mcs.anl.gov>
> wrote:
> >
> >    This behavior where the residual norm jumps at restart indicates
> something is very very wrong. Run with the option
> -ksp_monitor_true_residual and I think you'll see the true residual is not
> decreasing as is the preconditioned residual. My guess is that your "action
> of the matrix" is incorrect and not actually a linear operator. Try using
> MatComputeExplicitOperator() and see what explicit matrix it produces, is
> it what you expect?
> >
> >    Barry
> >
> >
> >
> >
> > > On May 27, 2019, at 11:33 PM, Jan Izak Cornelius Vermaak via
> petsc-users <petsc-users at mcs.anl.gov> wrote:
> > >
> > > Hi all,
> > >
> > > So I am faced with this debacle. I have a neutron transport solver
> with a sweep code that can compute the action of the matrix on a vector.
> > >
> > > I use a matrix shell to set up the action of the matrix. The method
> works but only if I can get the solution converged before GMRES restarts.
> It gets the right answer. Now my first problem is (and I only saw this when
> I hit the first restart) is that it looks like the solver completely resets
> after the GMRES-restart. Below is an iteration log with restart interval
> set to 10. At first I thought it wasn't updating the initial guess but it
> became clear that it initial guess always had no effect. I do set
> KSPSetInitialGuessNonZero but it has no effect.
> > >
> > > Is the matrix-free business defaulting my initial guess to zero
> everytime? What can I do to actually supply an initial guess? I've used
> PETSc for diffusion many times and the initial guess always works, just not
> now.
> > >
> > > [0]  Computing b
> > > [0]  Iteration 0 Residual 169.302
> > > [0]  Iteration 1 Residual 47.582
> > > [0]  Iteration 2 Residual 13.2614
> > > [0]  Iteration 3 Residual 4.46795
> > > [0]  Iteration 4 Residual 1.03038
> > > [0]  Iteration 5 Residual 0.246807
> > > [0]  Iteration 6 Residual 0.0828341
> > > [0]  Iteration 7 Residual 0.0410627
> > > [0]  Iteration 8 Residual 0.0243749
> > > [0]  Iteration 9 Residual 0.0136067
> > > [0]  Iteration 10 Residual 169.302
> > > [0]  Iteration 11 Residual 47.582
> > > [0]  Iteration 12 Residual 13.2614
> > > [0]  Iteration 13 Residual 4.46795
> > > [0]  Iteration 14 Residual 1.03038
> > > [0]  Iteration 15 Residual 0.246807
> > > [0]  Iteration 16 Residual 0.0828341
> > > [0]  Iteration 17 Residual 0.0410627
> > > [0]  Iteration 18 Residual 0.0243749
> > > [0]  Iteration 19 Residual 0.0136067
> > > [0]  Iteration 20 Residual 169.302
> > >
> > > --
> > > Jan Izak Cornelius Vermaak
> > > (M.Eng Nuclear)
> > > Email: janicvermaak at gmail.com
> > > Cell:    +1-979-739-0789
> >
> >
> >
> > --
> > Jan Izak Cornelius Vermaak
> > (M.Eng Nuclear)
> > Email: janicvermaak at gmail.com
> > Cell:    +1-979-739-0789
>
>

-- 
Jan Izak Cornelius Vermaak
(M.Eng Nuclear)
Email: janicvermaak at gmail.com
Cell:    +1-979-739-0789
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20190529/44891acb/attachment-0001.html>