[petsc-dev] block size information in gamg

[Barry Smith]

  Stefano,

   There is currently a PCFieldSplitSetBlockSize() and a PCSPAISetBlockSize()

PetscErrorCode MatSetBlockSize(Mat mat,PetscInt bs)
{
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
  PetscValidLogicalCollectiveInt(mat,bs,2);
  ierr = PetscLayoutSetBlockSize(mat->rmap,bs);CHKERRQ(ierr);
  ierr = PetscLayoutSetBlockSize(mat->cmap,bs);CHKERRQ(ierr);

PetscErrorCode PetscLayoutSetBlockSize(PetscLayout map,PetscInt bs)
{
  PetscFunctionBegin;
  if (bs < 0) PetscFunctionReturn(0);
  if (map->n > 0 && map->n % bs) SETERRQ2(map->comm,PETSC_ERR_ARG_INCOMP,"Local size %D not compatible with block size %D",map->n,bs);
  if (map->range && map->bs > 0 && map->bs != bs) SETERRQ2(map->comm,PETSC_ERR_ARG_INCOMP,"Cannot change block size %D to %D",map->bs,bs);
  if (map->mapping) {
    PetscInt       lbs;
    PetscErrorCode ierr;

    ierr = ISLocalToGlobalMappingGetBlockSize(map->mapping,&lbs);CHKERRQ(ierr);
    if (lbs != bs) SETERRQ2(map->comm,PETSC_ERR_PLIB,"Blocksize of localtoglobalmapping %D must match that of layout %D",lbs,bs);
  }
  map->bs = bs;
  PetscFunctionReturn(0);

So it doesn't want to change the block size if map->range has been set (i.e. PetscLayoutSetUp()) was called. But if the matrix is AIJ so long as  !(map->n % bs) one can actually set the block size at any time. So how about if I change the Mat code to allow setting the block size late for AIJ?  Maybe we can also get rid of PCFieldSplitSetBlockSize() and a PCSPAISetBlockSize()?

One issue that comes up with this change is: should setting the block size later automatically reset the block size of the two inner AIJ matrices? Or should they remain one? I could have it update the inner block sizes at the expense of a tiny bit more code. But one complication is that MatSetUpMultiply_MPIAIJ() always sets the column block size of the off-diagaonl matrix to 1; so this needs to be respected.

The PETSc design policy of "there is one way to do something" dictates that it is best not to have PCSetBlockSize() unless absolutely necessary so I am trying to see if it is absolutely necessary.

  Barry




> On Oct 10, 2016, at 10:18 AM, Stefano Zampini <stefano.zampini at gmail.com> wrote:
> 
> Mark,
> 
> I was wondering if you could provide an API/command line customization for PCGAMG to explicitly provide the "block size" of the problem to the gamg solver, instead of using the block size of the matrix. The matrix block size could be always used as a fallback if this new information as not been provided.
> 
> I'm interfacing a complicated finite element package to PETSc , where you cannot do any assumption on the blocksize of the problem during matrix setup, and (AFAIK) I cannot change the block size with a call to MatSetBlockSize after the matrix has been finalized without hacking (that I would like to avoid). 
> 
> Thanks
> -- 
> Stefano