[petsc-users] Question about MatMultDiagonalBlock

[Barry Smith]

On Oct 5, 2012, at 7:38 PM, John Fettig <john.fettig at gmail.com> wrote:

> In mpiaij.c, there is this if statement:
> 
>   if (a->inode.size) mat->ops->multdiagonalblock = MatMultDiagonalBlock_MPIAIJ;
> 
> which prevents the operation MatMultDiagonalBlock from being set if
> a->inode.size is zero.  Why is this?  What is wrong with the diagonal
> block when inode.size==0?
> 
> Thanks,
> John

   John,

     When innode.size is zero the "diagonal blocks" are all 1 by 1.  So, yes, there is really no theoretical reason we couldn't provide a function for that case. The current MatMultDiagonalBlock_SeqAIJ_Inode() uses the inode data structure to do the computation so could not be used directly in that case. 

     MatMultDiagonalBlock is only used for Eisenstat's trick for AIJ matrices with inodes and when it has no inodes it used VecPointWiseMult() with the diagonal of the matrix.

#undef __FUNCT__
#define __FUNCT__ "PCApply_Eisenstat"
static PetscErrorCode PCApply_Eisenstat(PC pc,Vec x,Vec y)
{
  PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;
  PetscErrorCode ierr;
  PetscBool      hasop;

  PetscFunctionBegin;
  if (eis->usediag)  {
    ierr = MatHasOperation(pc->pmat,MATOP_MULT_DIAGONAL_BLOCK,&hasop);CHKERRQ(ierr);
    if (hasop) {
      ierr = MatMultDiagonalBlock(pc->pmat,x,y);CHKERRQ(ierr);
    } else {
      ierr = VecPointwiseMult(y,x,eis->diag);CHKERRQ(ierr);
    }
  } else {ierr = VecCopy(x,y);CHKERRQ(ierr);}
  PetscFunctionReturn(0);
}

  Are you using MatMultDiagonalBlock() for some other case?

   Barry