<div dir="ltr">On Mon, Apr 15, 2013 at 4:08 PM, Hugo Gagnon <span dir="ltr"><<a href="mailto:opensource.petsc@user.fastmail.fm" target="_blank">opensource.petsc@user.fastmail.fm</a>></span> wrote:<br><div class="gmail_extra">
<div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div style="word-wrap:break-word">For the problem I'm describing my serial in-house solver does not work with ILU(0) but works with ILU(3). I have no option to run Jacobi. When I apply the same problem to PETSc's PC solver with ILU(3) in serial I get KSP_DIVERGED_INDEFINITE_PC on the first iteration (in MPI the solution somewhat converges but very slowly).</div>
</blockquote><div><br></div><div style>As Mark said, ILU(3) does not preserve either symmetry or definiteness.</div><div style><br></div><div style> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div style="word-wrap:break-word"><div>call KSPGetPC(Pksp,Ppc,Pierr)</div><div>call PCSetType(Ppc,PCILU,Pierr)</div><div>call PCFactorSetLevels(Ppc,3,Pierr)</div><div><br></div><div>This effectively changes the fill level from 0 to 3, right?</div>
<div><br></div><div><div><div>
<div style="text-indent:0px;letter-spacing:normal;font-variant:normal;text-align:-webkit-auto;font-style:normal;font-weight:normal;line-height:normal;text-transform:none;font-size:medium;white-space:normal;font-family:Helvetica;word-wrap:break-word;word-spacing:0px">
<div>--</div><div> Hugo Gagnon</div></div>
</div>
<br><div><div>On 2013-04-15, at 4:35 PM, Mark F. Adams <<a href="mailto:mark.adams@columbia.edu" target="_blank">mark.adams@columbia.edu</a>> wrote:</div><br><blockquote type="cite">ILU is not guaranteed to stay positive and CG requires this. PETSc's ILU is ILU(0). If your in-house PCG solver works with ILU(0) then there is probably a bug in your construction of the PETSc matrix. If your in-house code does not use ILU(0) then I would try PCJACOBI and verify that you get the same residaul history as your in-house solver (assuming that you can do Jacobi).<br>
<br><br>On Apr 15, 2013, at 4:04 PM, Jed Brown <<a href="mailto:jedbrown@mcs.anl.gov" target="_blank">jedbrown@mcs.anl.gov</a>> wrote:<br><br><blockquote type="cite">Hugo Gagnon <<a href="mailto:opensource.petsc@user.fastmail.fm" target="_blank">opensource.petsc@user.fastmail.fm</a>> writes:<br>
<br><blockquote type="cite">Hi,<br><br>I hope you don't mind me writing to you off list, I just don't think<br>my problem would be helpful to others anyway.<br></blockquote><br>I'd rather keep discussions on the list. We (PETSc developers) cannot<br>
scale to have private conversations with all users.<br><br><blockquote type="cite">I'm trying to converge a problem with KSPCG, which I know for sure<br>that both A and b, as in Ax=b, are correct. I can successfully<br>
converge this particular problem using our own serial in-house PCG<br>solver, whereas with PETSc the solution blows up at the first<br>iteration with error code -8.<br></blockquote><br>-ksp_converged_reason prints why, or use KSPConvergedReasons to get the string.<br>
<br> KSP_DIVERGED_INDEFINITE_PC = -8,<br><br><blockquote type="cite"><br>I've included parts of my code below. Again, I've triple checked that<br>I build A and b correctly. Being a beginner with PETSc and parallel<br>
solvers in general, I'd appreciate if you could give me some simple<br>pointers as to why my problem is blowing up and how can I improve my<br>code.<br></blockquote><br><a href="http://scicomp.stackexchange.com/questions/513/why-is-my-iterative-linear-solver-not-converging" target="_blank">http://scicomp.stackexchange.com/questions/513/why-is-my-iterative-linear-solver-not-converging</a><br>
<br><blockquote type="cite">Note that this code does work for easier problems, although I've<br>noticed that in general PETSc takes way more iterations to converge<br>than our in-house solver. Would you have an idea as to why this is<br>
so?<br></blockquote><br>You are either doing a different algorithm or evaluating convergence<br>differently. What preconditioner are you using in serial? Why aren't<br>you using PETSc's multigrid? Send the diagnostics described in the link<br>
if you want more help.<br><br><blockquote type="cite">I know I'm asking quite a bit but I've been struggling with this<br>problem for what I feel far too long!<br><br>Thank you for your help,<br>--<br> Hugo Gagnon<br>
<br> !-- initialize petsc with the current communicator<br> PETSC_COMM_WORLD = COMM_CURRENT<br> call PetscInitialize(PETSC_NULL_CHARACTER,Pierr)<br><br> !-- instantiate the lhs and solution vectors<br> call VecCreate(PETSC_COMM_WORLD,Pb,Pierr)<br>
call VecSetSizes(Pb,PETSC_DECIDE,numVar,Pierr))<br> call VecSetType(Pb,VECSTANDARD,Pierr)<br> call VecDuplicate(Pb,Psol,Pierr)<br><br> !-- instantiate the rhs<br> call MatCreate(PETSC_COMM_WORLD,Paa,Pierr)<br> call MatSetSizes(Paa,PETSC_DECIDE,PETSC_DECIDE,numVar,numVar,Pierr)<br>
call MatSetType(Paa,MATAIJ,Pierr) !-- csr format<br><br> !-- preallocate memory<br> mnnz = 81 !-- knob<br> call MatSeqAIJSetPreallocation(Paa,mnnz,PETSC_NULL_INTEGER,Pierr)<br> call MatMPIAIJSetPreallocation(Paa,mnnz,PETSC_NULL_INTEGER,mnnz,PETSC_NULL_INTEGER,Pierr)<br>
<br> !-- instantiate the linear solver context<br> call KSPCreate(PETSC_COMM_WORLD,Pksp,Pierr)<br> call KSPSetType(Pksp,KSPCG,Pierr)<br><br> !-- use a custom monitor function<br> call KSPMonitorSet(Pksp,KSPMonitorPETSc,PETSC_NULL_OBJECT,PETSC_NULL_FUNCTION,Pierr)<br>
<br> !-- instantiate a scatterer so that root has access to full Psol<br> call VecScatterCreateToZero(Psol,Psols,Psol0,Pierr)<br><br> call MatZeroEntries(Paa,Pierr)<br><br> *** build Psol, Pb and Paa ****<br><br> !-- assemble (distribute) the initial solution, the lhs and the rhs<br>
call VecAssemblyBegin(Psol,Pierr)<br> call VecAssemblyBegin(Pb,Pierr)<br> call MatAssemblyBegin(Paa,MAT_FINAL_ASSEMBLY,Pierr)<br> call VecAssemblyEnd(Psol,Pierr)<br> call VecAssemblyEnd(Pb,Pierr)<br> call MatAssemblyEnd(Paa,MAT_FINAL_ASSEMBLY,Pierr)<br>
call MatSetOption(Paa,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE,Pierr)<br><br> !-- setup pcg, reusing the lhs as preconditioner<br> call KSPSetOperators(Pksp,Paa,Paa,SAME_NONZERO_PATTERN,Pierr)<br> call KSPSetTolerances(Pksp,tol,PETSC_DEFAULT_DOUBLE_PRECISION,PETSC_DEFAULT_DOUBLE_PRECISION,grid%pcg_its,Pierr)<br>
<br> !-- setup the preconditioner with the desired fill level<br> call KSPGetPC(Pksp,Ppc,Pierr)<br> if (nProc > 1) then<br> call PCSetType(Ppc,PCBJACOBI,Pierr)<br> call KSPSetup(Pksp,Pierr)<br> call PCBJacobiGetSubKSP(Ppc,j,j,Psubksp,Pierr)<br>
call KSPGetPC(Psubksp(1),Ppc,Pierr)<br> end if<br> call PCSetType(Ppc,PCILU,Pierr)<br> call PCFactorSetLevels(Ppc,grid%pcg_fill,Pierr)<br><br> !-- do not zero out the solution vector<br> call KSPSetInitialGuessNonzero(Pksp,PETSC_TRUE,Pierr)<br>
call KSPSetFromOptions(Pksp,Pierr)<br><br> !-- let petsc do its magic<br> call KSPSolve(Pksp,Pb,Psol,Pierr)<br> call KSPGetConvergedReason(Pksp,Pconv,Pierr)<br></blockquote><br></blockquote><br></blockquote></div>
<br></div></div></div></blockquote></div><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
</div></div>