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