[petsc-users] KSPCG solution blow up

[Mark F. Adams]

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).


On Apr 15, 2013, at 4:04 PM, Jed Brown <jedbrown at mcs.anl.gov> wrote:

> Hugo Gagnon <opensource.petsc at user.fastmail.fm> writes:
> 
>> Hi,
>> 
>> I hope you don't mind me writing to you off list, I just don't think
>> my problem would be helpful to others anyway.
> 
> I'd rather keep discussions on the list.  We (PETSc developers) cannot
> scale to have private conversations with all users.
> 
>> I'm trying to converge a problem with KSPCG, which I know for sure
>> that both A and b, as in Ax=b, are correct.  I can successfully
>> converge this particular problem using our own serial in-house PCG
>> solver, whereas with PETSc the solution blows up at the first
>> iteration with error code -8.
> 
> -ksp_converged_reason prints why, or use KSPConvergedReasons to get the string.
> 
>               KSP_DIVERGED_INDEFINITE_PC       = -8,
> 
>> 
>> I've included parts of my code below.  Again, I've triple checked that
>> I build A and b correctly.  Being a beginner with PETSc and parallel
>> solvers in general, I'd appreciate if you could give me some simple
>> pointers as to why my problem is blowing up and how can I improve my
>> code.
> 
> http://scicomp.stackexchange.com/questions/513/why-is-my-iterative-linear-solver-not-converging
> 
>> Note that this code does work for easier problems, although I've
>> noticed that in general PETSc takes way more iterations to converge
>> than our in-house solver.  Would you have an idea as to why this is
>> so?
> 
> You are either doing a different algorithm or evaluating convergence
> differently.  What preconditioner are you using in serial?  Why aren't
> you using PETSc's multigrid?  Send the diagnostics described in the link
> if you want more help.
> 
>> I know I'm asking quite a bit but I've been struggling with this
>> problem for what I feel far too long!
>> 
>> Thank you for your help,
>> --
>>  Hugo Gagnon
>> 
>>    !-- initialize petsc with the current communicator
>>    PETSC_COMM_WORLD = COMM_CURRENT
>>    call PetscInitialize(PETSC_NULL_CHARACTER,Pierr)
>> 
>>    !-- instantiate the lhs and solution vectors
>>    call VecCreate(PETSC_COMM_WORLD,Pb,Pierr)
>>    call VecSetSizes(Pb,PETSC_DECIDE,numVar,Pierr))
>>    call VecSetType(Pb,VECSTANDARD,Pierr)
>>    call VecDuplicate(Pb,Psol,Pierr)
>> 
>>    !-- instantiate the rhs
>>    call MatCreate(PETSC_COMM_WORLD,Paa,Pierr)
>>    call MatSetSizes(Paa,PETSC_DECIDE,PETSC_DECIDE,numVar,numVar,Pierr)
>>    call MatSetType(Paa,MATAIJ,Pierr) !-- csr format
>> 
>>    !-- preallocate memory
>>    mnnz = 81 !-- knob
>>    call MatSeqAIJSetPreallocation(Paa,mnnz,PETSC_NULL_INTEGER,Pierr)
>>    call MatMPIAIJSetPreallocation(Paa,mnnz,PETSC_NULL_INTEGER,mnnz,PETSC_NULL_INTEGER,Pierr)
>> 
>>    !-- instantiate the linear solver context
>>    call KSPCreate(PETSC_COMM_WORLD,Pksp,Pierr)
>>    call KSPSetType(Pksp,KSPCG,Pierr)
>> 
>>    !-- use a custom monitor function
>>    call KSPMonitorSet(Pksp,KSPMonitorPETSc,PETSC_NULL_OBJECT,PETSC_NULL_FUNCTION,Pierr)
>> 
>>    !-- instantiate a scatterer so that root has access to full Psol
>>    call VecScatterCreateToZero(Psol,Psols,Psol0,Pierr)
>> 
>>    call MatZeroEntries(Paa,Pierr)
>> 
>>    *** build Psol, Pb and Paa ****
>> 
>>    !-- assemble (distribute) the initial solution, the lhs and the rhs
>>    call VecAssemblyBegin(Psol,Pierr)
>>    call VecAssemblyBegin(Pb,Pierr)
>>    call MatAssemblyBegin(Paa,MAT_FINAL_ASSEMBLY,Pierr)
>>    call VecAssemblyEnd(Psol,Pierr)
>>    call VecAssemblyEnd(Pb,Pierr)
>>    call MatAssemblyEnd(Paa,MAT_FINAL_ASSEMBLY,Pierr)
>>    call MatSetOption(Paa,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE,Pierr)
>> 
>>    !-- setup pcg, reusing the lhs as preconditioner
>>    call KSPSetOperators(Pksp,Paa,Paa,SAME_NONZERO_PATTERN,Pierr)
>>    call KSPSetTolerances(Pksp,tol,PETSC_DEFAULT_DOUBLE_PRECISION,PETSC_DEFAULT_DOUBLE_PRECISION,grid%pcg_its,Pierr)
>> 
>>    !-- setup the preconditioner with the desired fill level
>>    call KSPGetPC(Pksp,Ppc,Pierr)
>>    if (nProc > 1) then
>>      call PCSetType(Ppc,PCBJACOBI,Pierr)
>>      call KSPSetup(Pksp,Pierr)
>>      call PCBJacobiGetSubKSP(Ppc,j,j,Psubksp,Pierr)
>>      call KSPGetPC(Psubksp(1),Ppc,Pierr)
>>    end if
>>    call PCSetType(Ppc,PCILU,Pierr)
>>    call PCFactorSetLevels(Ppc,grid%pcg_fill,Pierr)
>> 
>>    !-- do not zero out the solution vector
>>    call KSPSetInitialGuessNonzero(Pksp,PETSC_TRUE,Pierr)
>>    call KSPSetFromOptions(Pksp,Pierr)
>> 
>>    !-- let petsc do its magic
>>    call KSPSolve(Pksp,Pb,Psol,Pierr)
>>    call KSPGetConvergedReason(Pksp,Pconv,Pierr)
>