[petsc-users] understanding gmres output

Tue May 1 07:13:22 CDT 2018

On Tue, May 1, 2018 at 7:15 AM, Ganesh Diwan <gcdiwan at gmail.com> wrote:

> You may/likely won't get the same convergence with Matlab not due to GMRES
>> but due to differences in the ILU in Matlab and PETSc. Better to run with
>> Jacobi or no preconditioner if you wish to get consistent results.
>
>
> Not sure how to specify 'no' preconditioner
>

-pc_type none

> but I tried ksp.setOperators(A,Id) where Id is the nxn identity matrix and
> this gave me consistent results with Matlab.
>
> I am working on the Helmholtz problems with finite elements
>

Notoriously difficult. For example
https://www.unige.ch/~gander/Preprints/HelmholtzReview.pdf

> and I use a preconditioning matrix which is supposed to help with the
> GMRES convergence. From the options set in the code in my previous post, I
> believe PETSc will compute the norm of preconditioned residual (default for
> left preconditioning) and use it as the stopping criterion. From the
> Matlab documentation, I understand that it will use the norm of the
> preconditioned residual (same as PETSc). Despite using the same
> parameters (number of restarts, maximum iterations, tolerance and the same
> preconditioning matrix), I see PETSc and Matlab diverge in terms of the
> GMRES iterations as I refine my finite element mesh (I get consistent
> results only up to a certain mesh refinement level). Could you comment?
>

The problem becomes more ill-conditioned as you refine, and thus more
sensitive to small numerical errors such
as round-off. This sensitivity results in different answers when equivalent
computations are done in a different order.
Moreover, its possible that PETSc and Matlab for using different
orthogonalization routines. Classical Gram-Schmidt,
modified Gram-Schmidt, and Householder will likely give different answers.

  Thanks,

     Matt

>
>  Thank you, Ganesh
>
> On Thu, Apr 26, 2018 at 7:24 PM, Smith, Barry F. <bsmith at mcs.anl.gov>
> wrote:
>
>>
>> *  You may/likely won't get the same convergence with Matlab not due to
>> GMRES but due to differences in the ILU in Matlab and PETSc. Better to run
>> with Jacobi or no preconditioner if you wish to get consistent results.
>>
>> *   Check carefully if Matlab using left or right preconditioning for
>> GMRES (PETSc uses left by default).
>>
>> > I see that the residual < norm(b)*rtol for the 4th iteration itself, I
>> am not sure then why GMRES continues for one more iteration.
>>
>> * Since PETSc is using left preconditioning it is norm(B*b)*rtol that
>> determines the convergence, not norm(b)*rtol. Where B is application of
>> preconditioner.
>>
>> * The residual history for GMRES includes everything through any restarts
>> so if it takes 90 total iterations (and the restart is 30) the residual
>> history should have 90 values.
>>
>> > In the above specific example, it took GMRES a total of 5 iterations to
>> converge. Does it mean the convergence was achieved without having to
>> restart?
>>
>>   In most cases if you need to restart then you need a better
>> preconditioner.
>>
>>
>>    Barry
>>
>>
>>
>>
>> > On Apr 26, 2018, at 9:46 AM, Ganesh Diwan <gcdiwan at gmail.com> wrote:
>> >
>> > Dear Petsc developers,
>> >
>> > I am trying to learn to use PETSc GMRES with petsc4py and also trying
>> to follow the codes in petsc/src/ksp/ksp/examples/tutorials. Here is the
>> Python code I use to read a linear system and solve it with GMRES in PETSc
>> (petsc4py)
>> >
>> > # gmrestest.py
>> > # usual imports
>> > import numpy as np
>> > import scipy.io as sio
>> > import scipy
>> > from sys import getrefcount
>> > from petsc4py import PETSc
>> >
>> > # read the matrix data
>> > mat_contents = sio.loadmat('data.mat')
>> > mat_A = mat_contents['A']
>> > vec_b = mat_contents['b']
>> > n = mat_contents['n']
>> > vec_x = mat_contents['x']
>> > mat_nnzA = mat_contents['nnzA']
>> >
>> > # form the petsc matrices
>> > x = PETSc.Vec().createWithArray(vec_x)
>> > b = PETSc.Vec().createWithArray(vec_b)
>> > p1=mat_A.indptr
>> > p2=mat_A.indices
>> > p3=mat_A.data
>> > A = PETSc.Mat().createAIJ(size=mat_A.shape,csr=(p1,p2,p3))
>> > A = scipy.transpose(A) # transpose the csr format as my matrix is
>> column major originally
>> > A.assemblyBegin()
>> > A.assemblyEnd()
>> >
>> > # solve
>> > ksp = PETSc.KSP()
>> > ksp.create(PETSc.COMM_WORLD)
>> > ksp.setType('gmres')
>> > ksp.setOperators(A)
>> > ksp.setFromOptions()
>> > rtol = 1e-4
>> > ksp.setTolerances(rtol, 1e-5, 10000., 10000)
>> > ksp.view()
>> > ksp.setConvergenceHistory()
>> > ksp.solve(b, x)
>> >
>> > # print
>> > print 'iterations = ', ksp.getIterationNumber()
>> > print 'residual = ', '{:.2e}'.format(ksp.getResidualNorm())#  %.2E
>> ksp.getResidualNorm()
>> > print 'norm(b)*rtol =  ', '{:.2e}'.format(PETSc.Vec.norm(b)*rtol)#
>> %.2E PETSc.Vec.norm(b)*rtol
>> > print 'converge reason# = ', ksp.getConvergedReason()
>> > print 'residuals at each iter = ', ksp.getConvergenceHistory()
>> > #
>> >
>> > Here is the output from the above code for a linear system of dimension
>> 100 by 100.
>> >
>> >
>> > KSP Object: 1 MPI processes
>> >   type: gmres
>> >     restart=30, using Classical (unmodified) Gram-Schmidt
>> Orthogonalization with no iterative refinement
>> >     happy breakdown tolerance 1e-30
>> >   maximum iterations=10000, initial guess is zero
>> >   tolerances:  relative=0.0001, absolute=1e-05, divergence=10000.
>> >   left preconditioning
>> >   using DEFAULT norm type for convergence test
>> > PC Object: 1 MPI processes
>> >   type: ilu
>> >   PC has not been set up so information may be incomplete
>> >     out-of-place factorization
>> >     0 levels of fill
>> >     tolerance for zero pivot 2.22045e-14
>> >     matrix ordering: natural
>> >   linear system matrix = precond matrix:
>> >   Mat Object: 1 MPI processes
>> >     type: seqaij
>> >     rows=100, cols=100
>> >     total: nonzeros=2704, allocated nonzeros=2704
>> >     total number of mallocs used during MatSetValues calls =0
>> >       using I-node routines: found 25 nodes, limit used is 5
>> > iterations =  5
>> > residual =  1.38e-03
>> > norm(b)*rtol =   1.99e-01
>> > converge reason# =  2
>> > residuals at each iter =  [  2.05677686e+01   4.97916031e+00
>>  4.82888782e-01   1.16849581e-01
>> >    8.87159777e-03   1.37992327e-03]
>> >
>> > Sorry if this sounds stupid, but I am trying to understand the output
>> from PETSc by contrasting it with Matlab GMRES.
>> > I see that the residual < norm(b)*rtol for the 4th iteration itself, I
>> am not sure then why GMRES continues for one more iteration. Secondly, how
>> does one get total number of iterations? For eg. let's say if it takes 3
>> outer iterations each with the default restart of 30, then one would expect
>> the length of residual vector to be 150. In the above specific example, it
>> took GMRES a total of 5 iterations to converge. Does it mean the
>> convergence was achieved without having to restart?
>> >
>> > Thank you, Ganesh
>>
>>
>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ <http://www.caam.rice.edu/~mk51/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20180501/ecf4a824/attachment-0001.html>