Il 31/08/21 17:32, Jed Brown ha scritto:
Matteo Semplice <[email protected]> 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
