# Matrix free example snes/ex20.c

Barry Smith bsmith at mcs.anl.gov
Sun Dec 30 16:51:20 CST 2007

```   Vijay,

This is a very cool problem.

Because of the exact symmetry of the domain the EXACT Jacobian
at each step has exactly 9 different eigenvalues. This means the
GMRES will take exactly 9 iterations (and "completely" converge in
the ninth iteration) if the "exact" Jacobian is used.  You can run with
-pc_type none -snes_mf -ksp_monitor_singular_value -
ksp_plot_eigenvalues -display :0.0 -draw_pause -1
to see the 9 eigenvalues.

Now run without the -snes_mf option. You will see the first Newton
iteration's
eigenvalues still look like 9; but starting at the second Newton
iteration the
"identical" eigenvalues are now not all identically placed so GMRES
needs
more iterations.

The question then becomes how come the matrix-free application of
the Jacobian is more accurate than actually computing it as a sparse
matrix then applying it? Here is my non-rigorous answer; the
multiplication
of the sparse matrix values (even if very accurate) against the vector
introduces
some rounding error that screws up the eigenvalues slightly. For some
reason
for this problem the matrix-free application is accurate enough not to
perturb the eigenvalues.

Barry

On Dec 30, 2007, at 2:19 PM, Vijay M wrote:

> I ran both the cases with -ksp_monitor on and have attached the
> output in
> two different files. 1.txt is the Jfree case and 2.txt is the
> analytical
> case.
>
> Barry, the example problem is ex20 from the snes tutorial directory.
> The
> petsc version is 2.3.3-p7 if that helps to clear things a little.
> Now I
> haven't yet completely checked for a bug in the analytical Jacobian
> but I
> would imagine that if it were incorrect, wouldn't that affect only
> how the
> nonlinear iteration converges and not the linear iteration since the
> matrix
> sparsity structure is still the same (well assuming the condition
> number is
> not very different from the exact Jacobian !). Just my 2 cents.
>
> Anyway, I will look into the code for ex20 and then see if something
> is
> messed up. Let me know if you find out the problem from the output.
>
> Thanks,
> Vijay
>
>> I understand that both the methods will not give me the same number
>> of total
>> linear iterations but a factor of 2 seems a little odd to me.
>
>    Yes, this is surprising.
>
>    Run with -ksp_monitor how are the linear convergence different?
>
>> another question whether the user can actually change the epsilon
>> used for
>> computing the perturbation in J-free scheme or is this fixed in
>> PETSc ?
>
>    Yes, see the manual page for MatMFFDSetFromOptions() and related
> manual
> pages.
>
>>
>>
>> If not, then what do you think is the reason for this ?
>
>    Bug in your analytic Jacobian? Run with -snes_monitor and -
> ksp_monitor and
> send all output.
>
>    Barry
>
>> Do let me know your
>> comments when you get some time. Thanks.
>>
>> Vijay
>>
>> -----Original Message-----
>> From: owner-petsc-users at mcs.anl.gov [mailto:owner-petsc-users at mcs.anl.gov
>> ]
>> On Behalf Of Matthew Knepley
>> Sent: Saturday, December 29, 2007 9:05 PM
>> To: petsc-users at mcs.anl.gov
>> Subject: Re: Matrix free example snes/ex20.c
>>
>> On Dec 29, 2007 8:07 PM, Vijay M <vijay.m at gmail.com> wrote:
>>> Hi all,
>>>
>>> I was trying to compile and run the ex20.c example code in the
>>> tutorial
>>> section of SNES. Although it does not explicitly specify that -
>>> snes_mf
>>> option can be used, my understanding is that as long as a nonlinear
>> residual
>>> function is written correctly, PETSc will calculate via finite
>>> difference
>>> the action of the Jacobian on a given vector. Is that correct ?
>>
>> Yes.
>>
>>> Now if that is the case, then please observe the discrepancy in the
>>> number
>>> of linear iterations taken with an analytical Jacobian and matrix-
>>> free
>>> option. What puzzles me is that the SNES function norm are quite
>>> close for
>>> both the methods but the linear iterations differ by a factor of 3.
>>> Why
>>> exactly is this ?
>>
>> There is no PC when using -snes_mf whereas the default is ILU for the
>> analytic
>> Jacobian.
>>
>>  Matt
>>
>>> Here's the output to make this clearer.
>>>
>>> vijay :mpirun -np 1 ex20 -ksp_type gmres -snes_monitor
>>>
>>> 0 SNES Function norm 2.271442542876e-01
>>>
>>> 1 SNES Function norm 6.881516100891e-02
>>>
>>> 2 SNES Function norm 1.813939751552e-02
>>>
>>> 3 SNES Function norm 2.354176462207e-03
>>>
>>> 4 SNES Function norm 3.063728077362e-05
>>>
>>> 5 SNES Function norm 3.106106268946e-08
>>>
>>> 6 SNES Function norm 5.344742712545e-12
>>>
>>> 0 SNES Function norm 2.271442542876e-01
>>>
>>> 1 SNES Function norm 6.881516100891e-02
>>>
>>> 2 SNES Function norm 1.813939751552e-02
>>>
>>> 3 SNES Function norm 2.354176462207e-03
>>>
>>> 4 SNES Function norm 3.063728077362e-05
>>>
>>> 5 SNES Function norm 3.106106268946e-08
>>>
>>> 6 SNES Function norm 5.344742712545e-12
>>>
>>> Number of Newton iterations = 6
>>>
>>> Number of Linear iterations = 18
>>>
>>> Average Linear its / Newton = 3.000000e+00
>>>
>>>
>>>
>>> vijay :mpirun -np 1 ex20 -ksp_type gmres -snes_monitor -snes_mf
>>>
>>> 0 SNES Function norm 2.271442542876e-01
>>>
>>> 1 SNES Function norm 6.870629867542e-02
>>>
>>> 2 SNES Function norm 1.804335379848e-02
>>>
>>> 3 SNES Function norm 2.290074339682e-03
>>>
>>> 4 SNES Function norm 3.082384186373e-05
>>>
>>> 5 SNES Function norm 3.926396277038e-09
>>>
>>> 6 SNES Function norm 3.754922566585e-16
>>>
>>> 0 SNES Function norm 2.271442542876e-01
>>>
>>> 1 SNES Function norm 6.870629867542e-02
>>>
>>> 2 SNES Function norm 1.804335379848e-02
>>>
>>> 3 SNES Function norm 2.290074339682e-03
>>>
>>> 4 SNES Function norm 3.082384186373e-05
>>>
>>> 5 SNES Function norm 3.926396277038e-09
>>>
>>> 6 SNES Function norm 3.754922566585e-16
>>>
>>> Number of Newton iterations = 6
>>>
>>> Number of Linear iterations = 54
>>>
>>> Average Linear its / Newton = 9.000000e+00
>>>
>>>
>>>
>>> Thanks,
>>>
>>> Vijay
>>>
>>>
>>
>>
>>
>> --
>> What most experimenters take for granted before they begin their
>> experiments is infinitely more interesting than any results to which
>> -- Norbert Wiener
>>
> <1.txt><2.txt>

```