[petsc-users] Mat preallocation for SNES jacobian [WAS Re: Mat preallocation in case of variable stencils]

[Matteo Semplice]

Il 31/08/21 17:32, Jed Brown ha scritto:
> Matteo Semplice <matteo.semplice at uninsubria.it> writes:
>
>> Hi.
>>
>> We are writing a code for a FD scheme on an irregular domain and thus
>> the local stencil is quite variable: we have inner nodes, boundary nodes
>> and inactive nodes, each with their own stencil type and offset with
>> respect to the grid node. We currently create a matrix with
>> DMCreateMatrix on a DMDA and for now have set the option
>> MAT_NEW_NONZERO_LOCATIONS to PETSC_TRUE, but its time to render the code
>> memory-efficient. The layout created automatically is correct for inner
>> nodes, wrong for boundary ones (off-centered stencils) and redundant for
>> outer nodes.
>>
>> After the preprocessing stage (including stencil creation) we'd be in
>> position to set the nonzero pattern properly.
>>
>> Do we need to start from a Mat created by CreateMatrix? Or is it ok to
>> call DMCreateMatrix (so that the splitting among CPUs and the block size
>> are set by PETSc) and then call a MatSetPreallocation routine?
> You can call MatXAIJSetPreallocation after. It'll handle all matrix types so you don't have to shepherd data for all the specific preallocations.

Hi.

Actually I am still struggling with this... Let me explain.

My code relies on a SNES and the matrix I need to preallocate is the 
Jacobian. So I do:

in the main file
   ierr = DMCreateMatrix(ctx.daAll,&ctx.J);CHKERRQ(ierr);
   ierr = setJacobianPattern(ctx,ctx.J);CHKERRQ(ierr); //calling 
MatXAIJSetPreallocation on the second argument
   ierr = MatSetOption(ctx.J,MAT_NEW_NONZERO_LOCATIONS,*******); 
CHKERRQ(ierr);//allow new nonzeros?

   ierr = SNESSetFunction(snes,ctx.F      ,FormFunction,(void *) &ctx); 
CHKERRQ(ierr);
   ierr = SNESSetJacobian(snes,ctx.J,ctx.J,FormJacobian,(void *) &ctx); 
CHKERRQ(ierr);

   ierr = FormSulfationRHS(ctx, ctx.U0, ctx.RHS);CHKERRQ(ierr);
   ierr = SNESSolve(snes,ctx.RHS,ctx.U); CHKERRQ(ierr);

and

PetscErrorCode FormJacobian(SNES snes,Vec U,Mat J, Mat P,void *_ctx)

does (this is a 2 dof finite difference problem, so logically 2x2 blocks 
in J)

     ierr = setJac00(....,P) //calls to MatSetValues in the 00 block

     ierr = setJac01(....,P) //calls to MatSetValues in the 01 block

     ierr = setJac1X(....,P) //calls to MatSetValues in the 10 ans 11 block

     ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
     ierr = MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

If I run with MAT_NEW_NONZERO_LOCATIONS=TRUE, all runs fine and using 
the -info option I see that no mallocs are performed during Assembly.

Computing F
   0 SNES Function norm 7.672682917097e+02
Computing J
[0] MatAssemblyEnd_SeqAIJ(): Matrix size: 71874 X 71874; storage space: 
17661 unneeded,191714 used
[0] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 0
[0] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 27
[0] MatCheckCompressedRow(): Found the ratio (num_zerorows 
0)/(num_localrows 71874) < 0.6. Do not use CompressedRow routines.

If I omit the call to setJacobianPattern, info reports a nonsero number 
of mallocs, so somehow the setJacobianPattern routine should be doing 
its job correctly.

However, if I run with MAT_NEW_NONZERO_LOCATIONS=FALSE, the Jacobian is 
entirely zero and no error messages appear until the KSP tries to do its 
job:

Computing F
   0 SNES Function norm 7.672682917097e+02
Computing J
[0] MatAssemblyEnd_SeqAIJ(): Matrix size: 71874 X 71874; storage space: 
209375 unneeded,0 used
[0] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 0
[0] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 0
[0] MatCheckCompressedRow(): Found the ratio (num_zerorows 
71874)/(num_localrows 71874) > 0.6. Use CompressedRow routines.
... and then KSP complains!

I have tried adding MAT_FLUSH_ASSEMBLY calls inside the subroutines, but 
nothing changes.

So I have 2 questions:

1. If, as a temporary solution, I leave MAT_NEW_NONZERO_LOCATIONS=TRUE, 
am I going to incur in performance penalties even if no new nonzeros are 
created by my assembly routine?

2. Can you guess what is causing the problem?

Thanks

     Matteo