On Tue, May 1, 2018 at 7:15 AM, Ganesh Diwan <[email protected]> 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. <[email protected]> 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 <[email protected]> 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/>