<html><head><meta http-equiv="Content-Type" content="text/html charset=us-ascii"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; ">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><br></div><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 apple-content-edited="true">
<div style="color: rgb(0, 0, 0); font-family: Helvetica; font-size: medium; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; orphans: 2; text-align: -webkit-auto; text-indent: 0px; text-transform: none; white-space: normal; widows: 2; word-spacing: 0px; -webkit-text-size-adjust: auto; -webkit-text-stroke-width: 0px; word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; "><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">mark.adams@columbia.edu</a>> wrote:</div><br class="Apple-interchange-newline"><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">jedbrown@mcs.anl.gov</a>> wrote:<br><br><blockquote type="cite">Hugo Gagnon <<a href="mailto:opensource.petsc@user.fastmail.fm">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">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></body></html>