<div dir="ltr"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><span style="color:rgb(34,34,34);font-family:arial,sans-serif;font-size:12.8px;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">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.</span></blockquote><div> </div><div>Not sure how to specify 'no' preconditioner but I tried ksp.setOperators(A,Id) where Id is the nxn identity matrix and this gave me consistent results with Matlab.</div><div><br></div><div>I am working on the Helmholtz problems with finite elements 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. <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">From the Matlab documentation, I understand that it will use the norm of the preconditioned residual (same as PETSc). </span>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?</div><div><br></div><div> Thank you, Ganesh</div></div><div class="gmail_extra"><br><div class="gmail_quote">On Thu, Apr 26, 2018 at 7:24 PM, Smith, Barry F. <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><br>
* 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.<br>
<br>
* Check carefully if Matlab using left or right preconditioning for GMRES (PETSc uses left by default).<br>
<span class=""><br>
> 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.<br>
<br>
</span>* 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.<br>
<br>
* 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.<br>
<span class=""><br>
> 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?<br>
<br>
</span> In most cases if you need to restart then you need a better preconditioner.<br>
<span class="HOEnZb"><font color="#888888"><br>
<br>
Barry<br>
</font></span><div class="HOEnZb"><div class="h5"><br>
<br>
<br>
<br>
> On Apr 26, 2018, at 9:46 AM, Ganesh Diwan <<a href="mailto:gcdiwan@gmail.com">gcdiwan@gmail.com</a>> wrote:<br>
> <br>
> Dear Petsc developers,<br>
> <br>
> I am trying to learn to use PETSc GMRES with petsc4py and also trying to follow the codes in petsc/src/ksp/ksp/examples/<wbr>tutorials. Here is the Python code I use to read a linear system and solve it with GMRES in PETSc (petsc4py)<br>
> <br>
> # gmrestest.py<br>
> # usual imports<br>
> import numpy as np<br>
> import <a href="http://scipy.io" rel="noreferrer" target="_blank">scipy.io</a> as sio<br>
> import scipy<br>
> from sys import getrefcount<br>
> from petsc4py import PETSc<br>
> <br>
> # read the matrix data<br>
> mat_contents = sio.loadmat('data.mat')<br>
> mat_A = mat_contents['A']<br>
> vec_b = mat_contents['b']<br>
> n = mat_contents['n']<br>
> vec_x = mat_contents['x']<br>
> mat_nnzA = mat_contents['nnzA']<br>
> <br>
> # form the petsc matrices<br>
> x = PETSc.Vec().createWithArray(<wbr>vec_x)<br>
> b = PETSc.Vec().createWithArray(<wbr>vec_b)<br>
> p1=mat_A.indptr<br>
> p2=mat_A.indices<br>
> p3=mat_A.data<br>
> A = PETSc.Mat().createAIJ(size=<wbr>mat_A.shape,csr=(p1,p2,p3))<br>
> A = scipy.transpose(A) # transpose the csr format as my matrix is column major originally<br>
> A.assemblyBegin()<br>
> A.assemblyEnd()<br>
> <br>
> # solve<br>
> ksp = PETSc.KSP()<br>
> ksp.create(PETSc.COMM_WORLD)<br>
> ksp.setType('gmres')<br>
> ksp.setOperators(A)<br>
> ksp.setFromOptions()<br>
> rtol = 1e-4<br>
> ksp.setTolerances(rtol, 1e-5, 10000., 10000)<br>
> ksp.view()<br>
> ksp.setConvergenceHistory()<br>
> ksp.solve(b, x)<br>
> <br>
> # print<br>
> print 'iterations = ', ksp.getIterationNumber()<br>
> print 'residual = ', '{:.2e}'.format(ksp.<wbr>getResidualNorm())# %.2E ksp.getResidualNorm()<br>
> print 'norm(b)*rtol = ', '{:.2e}'.format(PETSc.Vec.<wbr>norm(b)*rtol)# %.2E PETSc.Vec.norm(b)*rtol<br>
> print 'converge reason# = ', ksp.getConvergedReason()<br>
> print 'residuals at each iter = ', ksp.getConvergenceHistory()<br>
> #<br>
> <br>
> Here is the output from the above code for a linear system of dimension 100 by 100.<br>
> <br>
> <br>
> KSP Object: 1 MPI processes<br>
> type: gmres<br>
> restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement<br>
> happy breakdown tolerance 1e-30<br>
> maximum iterations=10000, initial guess is zero<br>
> tolerances: relative=0.0001, absolute=1e-05, divergence=10000.<br>
> left preconditioning<br>
> using DEFAULT norm type for convergence test<br>
> PC Object: 1 MPI processes<br>
> type: ilu<br>
> PC has not been set up so information may be incomplete<br>
> out-of-place factorization<br>
> 0 levels of fill<br>
> tolerance for zero pivot 2.22045e-14<br>
> matrix ordering: natural<br>
> linear system matrix = precond matrix:<br>
> Mat Object: 1 MPI processes<br>
> type: seqaij<br>
> rows=100, cols=100<br>
> total: nonzeros=2704, allocated nonzeros=2704<br>
> total number of mallocs used during MatSetValues calls =0<br>
> using I-node routines: found 25 nodes, limit used is 5<br>
> iterations = 5<br>
> residual = 1.38e-03<br>
> norm(b)*rtol = 1.99e-01<br>
> converge reason# = 2<br>
> residuals at each iter = [ 2.05677686e+01 4.97916031e+00 4.82888782e-01 1.16849581e-01<br>
> 8.87159777e-03 1.37992327e-03]<br>
> <br>
> Sorry if this sounds stupid, but I am trying to understand the output from PETSc by contrasting it with Matlab GMRES.<br>
> 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?<br>
> <br>
> Thank you, Ganesh<br>
<br>
</div></div></blockquote></div><br></div>