<div dir="ltr">Dear Petsc developers,<div><br></div><div>I am trying to learn to use <span style="color:rgb(34,34,34);font-family:arial,sans-serif;font-size:small;font-style:normal;font-variant-ligatures:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;background-color:rgb(255,255,255);text-decoration-style:initial;text-decoration-color:initial;float:none;display:inline">PETSc</span> 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 <span style="color:rgb(34,34,34);font-family:arial,sans-serif;font-size:small;font-style:normal;font-variant-ligatures:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;background-color:rgb(255,255,255);text-decoration-style:initial;text-decoration-color:initial;float:none;display:inline">GMRES</span> in PETSc (petsc4py)</div><div><br></div><div># gmrestest.py</div><div><div># usual imports</div><div>import numpy as np</div><div>import <a href="http://scipy.io">scipy.io</a> as sio</div><div>import scipy</div><div>from sys import getrefcount</div><div>from petsc4py import PETSc</div><div><br></div><div># read the matrix data</div><div>mat_contents = sio.loadmat('data.mat')</div><div>mat_A = mat_contents['A']</div><div>vec_b = mat_contents['b']</div><div>n = mat_contents['n']</div><div>vec_x = mat_contents['x']</div><div>mat_nnzA = mat_contents['nnzA']</div><div><br></div><div># form the petsc matrices</div><div>x = PETSc.Vec().createWithArray(vec_x)</div><div>b = PETSc.Vec().createWithArray(vec_b)</div><div>p1=mat_A.indptr</div><div>p2=mat_A.indices</div><div>p3=mat_A.data</div><div>A = PETSc.Mat().createAIJ(size=mat_A.shape,csr=(p1,p2,p3))</div><div>A = scipy.transpose(A) # transpose the csr format as my matrix is column major originally</div><div>A.assemblyBegin()</div><div>A.assemblyEnd()</div><div><br></div><div># solve</div><div>ksp = PETSc.KSP()</div><div>ksp.create(PETSc.COMM_WORLD)</div><div>ksp.setType('gmres')</div><div>ksp.setOperators(A)</div><div>ksp.setFromOptions()</div><div>rtol = 1e-4</div><div>ksp.setTolerances(rtol, 1e-5, 10000., 10000)</div><div>ksp.view()</div><div>ksp.setConvergenceHistory()</div><div>ksp.solve(b, x)</div><div><br></div><div># print</div><div>print 'iterations = ', ksp.getIterationNumber()</div><div>print 'residual = ', '{:.2e}'.format(ksp.getResidualNorm())#  %.2E ksp.getResidualNorm()</div><div>print 'norm(b)*rtol =  ', '{:.2e}'.format(PETSc.Vec.norm(b)*rtol)# %.2E PETSc.Vec.norm(b)*rtol</div><div>print 'converge reason# = ', ksp.getConvergedReason()</div><div>print 'residuals at each iter = ', ksp.getConvergenceHistory()</div></div><div>#</div><div><br></div><div>Here is the output from the above code for a linear system of dimension 100 by 100.</div><div><br></div><div><br></div><div><div>KSP Object: 1 MPI processes<br></div><div>  type: gmres</div><div>    restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement</div><div>    happy breakdown tolerance 1e-30</div><div>  maximum iterations=10000, initial guess is zero</div><div>  tolerances:  relative=0.0001, absolute=1e-05, divergence=10000.</div><div>  left preconditioning</div><div>  using DEFAULT norm type for convergence test</div><div>PC Object: 1 MPI processes</div><div>  type: ilu</div><div>  PC has not been set up so information may be incomplete</div><div>    out-of-place factorization</div><div>    0 levels of fill</div><div>    tolerance for zero pivot 2.22045e-14</div><div>    matrix ordering: natural</div><div>  linear system matrix = precond matrix:</div><div>  Mat Object: 1 MPI processes</div><div>    type: seqaij</div><div>    rows=100, cols=100</div><div>    total: nonzeros=2704, allocated nonzeros=2704</div><div>    total number of mallocs used during MatSetValues calls =0</div><div>      using I-node routines: found 25 nodes, limit used is 5</div><div>iterations =  5</div><div>residual =  1.38e-03</div><div>norm(b)*rtol =   1.99e-01</div><div>converge reason# =  2</div><div>residuals at each iter =  [  2.05677686e+01   4.97916031e+00   4.82888782e-01   1.16849581e-01</div><div>   8.87159777e-03   1.37992327e-03]</div></div><div><br></div><div>Sorry if this sounds stupid, but I am trying to understand the output from PETSc by contrasting it with Matlab GMRES.</div><div>I see that the residual < norm(b)*rtol for the 4th iteration itself<span style="color:rgb(34,34,34);font-family:arial,sans-serif;font-size:small;font-style:normal;font-variant-ligatures:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;background-color:rgb(255,255,255);text-decoration-style:initial;text-decoration-color:initial;float:none;display:inline">, I am not sure then why GMRES continues for one more iteration. </span>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?</div><div><br></div><div>Thank you, Ganesh</div></div>