Forgive me if switching this to petsc-dev isn't appropriate, but I thought that I would mention that it is a (minor) possibility that this issue was caused by an instability in the way the Givens rotations are computed in PETSc's GMRES Hessenberg update (starting on line 390):<br>
<a href="http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/impls/gmres/gmres.c.html#KSPGMRES">http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/impls/gmres/gmres.c.html#KSPGMRES</a><br><br>Givens rotations are fairly tricky to avoid unnecessary over/underflow in, and it makes sense to use dlartg and friends, as they implement the algorithm described in this paper:<br>
<a href="http://www.netlib.org/lapack/lawns/lawn148.ps">http://www.netlib.org/lapack/lawns/lawn148.ps</a><br><br>Jack<br><br>On Sat, Jan 7, 2012 at 4:14 PM, Matthew Knepley <span dir="ltr"><<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>></span> wrote:<br>
<div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0pt 0pt 0pt 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div class="im">On Sat, Jan 7, 2012 at 4:00 PM, Mohamad M. Nasr-Azadani <span dir="ltr"><<a href="mailto:mmnasr@gmail.com" target="_blank">mmnasr@gmail.com</a>></span> wrote:<br>
</div><div class="gmail_quote"><div class="im"><blockquote class="gmail_quote" style="margin:0pt 0pt 0pt 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
Hi guys, <div><br></div><div>I am trying to narrow down an issue with my Poisson solver. </div><div>I have the following problem setup</div><div><br></div><div>Laplace(f) = rhs(x,z,y)</div><div>0 <= x,y,z <= (Lx,Ly,Lz)</div>



<div><br></div><div>I solve the Poisson equation in three dimensions with the analytical function f(x,y,z) defined by</div><div><br></div><div>f(x,z,y) = cos(2*pi*x/Lx)*cos(2*pi*y/Ly)*cos(2*pi*z/Lz) + K</div><div>where Lx = Ly =Lz = 1.0 and K is a constant I use to set f(Lx,Ly,Lz) = 0.0. </div>



<div><br></div><div>Second order descritization is used for the Poisson equation. </div><div>Also, Neumann boundary condition is used everywhere, but I set the top-right-front node's value to zero to get rid of the Nullspaced matrix manually. </div>



<div>I use 20 grid points in each direction. </div><div><br></div><div>The problem is: </div><div>I use GMRES(20) without any preconditioners (rtol = 1e-12) to solve the linear system. </div><div>It takes 77,000 iterations to converge!!!! </div>



<div><br></div><div>For the size of only 8,000 unknowns, even though the lsys is not preconditioned, I guess that is a LOT of iterations. </div><div>Next, I setup the exact same problem in MATLAB and use their GMRES solver function. </div>



<div>I set the same parameters and MATLAB tells me that it converges using only 3870 iterations. </div></blockquote><div><br></div></div><div>1) Matlab could be doing a lot of things. I am betting that they scale the problem, so -pc_type jacobi.</div>

<div><br></div><div>2) Why would anyone ever use GMRES without a preconditioner, particularly for a problem where several</div><div>    optimal PCs exist and are present in PETSc.</div><div><br></div><div>   Matt</div><div class="im">
<div>
 </div><blockquote class="gmail_quote" style="margin:0pt 0pt 0pt 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div>I know that there might be some internal differences between MATLAB and PETSc's implementations of this method, but given the fact that these two solvers are not preconditioned, I am wondering about this big difference? </div>



<div><br></div><div>Any ideas? </div><div><br></div><div>Best, </div><span><font color="#888888"><div>Mohamad</div><div><br></div>
</font></span></blockquote></div></div><span class="HOEnZb"><font color="#888888"><br><br clear="all"><div><br></div>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>

-- Norbert Wiener<br>
</font></span></blockquote></div><br>