[petsc-dev] what happened to RAP timers?

[Barry Smith]

#undef __FUNCT__
#define __FUNCT__ "MatPtAP"
/*@
   MatPtAP - Creates the matrix product C = P^T * A * P

   Neighbor-wise Collective on Mat

   Input Parameters:
+  A - the matrix
.  P - the projection matrix
.  scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX
-  fill - expected fill as ratio of nnz(C)/(nnz(A) + nnz(P)), use PETSC_DEFAULT if you do not have a good estimate
          if the result is a dense matrix this is irrelevent

   Output Parameters:
.  C - the product matrix

   Notes:
   C will be created and must be destroyed by the user with MatDestroy().

   This routine is currently only implemented for pairs of AIJ matrices and classes
   which inherit from AIJ.

   Level: intermediate

.seealso: MatPtAPSymbolic(), MatPtAPNumeric(), MatMatMult(), MatRARt()
@*/
PetscErrorCode MatPtAP(Mat A,Mat P,MatReuse scall,PetscReal fill,Mat *C)
{
  PetscErrorCode ierr;
  PetscErrorCode (*fA)(Mat,Mat,MatReuse,PetscReal,Mat*);
  PetscErrorCode (*fP)(Mat,Mat,MatReuse,PetscReal,Mat*);
  PetscErrorCode (*ptap)(Mat,Mat,MatReuse,PetscReal,Mat*)=NULL;
  PetscBool      viatranspose=PETSC_FALSE,viamatmatmatmult=PETSC_FALSE;

  PetscFunctionBegin;
  ierr = PetscOptionsGetBool(((PetscObject)A)->options,((PetscObject)A)->prefix,"-matptap_viatranspose",&viatranspose,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetBool(((PetscObject)A)->options,((PetscObject)A)->prefix,"-matptap_viamatmatmatmult",&viamatmatmatmult,NULL);CHKERRQ(ierr);

  PetscValidHeaderSpecific(A,MAT_CLASSID,1);
  PetscValidType(A,1);
  MatCheckPreallocated(A,1);
  if (!A->assembled) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
  if (A->factortype) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
  PetscValidHeaderSpecific(P,MAT_CLASSID,2);
  PetscValidType(P,2);
  MatCheckPreallocated(P,2);
  if (!P->assembled) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
  if (P->factortype) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

  if (A->rmap->N!= A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Matrix A must be square, %D != %D",A->rmap->N,A->cmap->N);
  if (P->rmap->N != A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",P->rmap->N,A->cmap->N);
  if (fill == PETSC_DEFAULT || fill == PETSC_DECIDE) fill = 2.0;
  if (fill < 1.0) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Expected fill=%g must be >= 1.0",(double)fill);

  if (scall == MAT_REUSE_MATRIX) {
    PetscValidPointer(*C,5);
    PetscValidHeaderSpecific(*C,MAT_CLASSID,5);
    if (viatranspose || viamatmatmatmult) {
      Mat Pt;
      ierr = MatTranspose(P,MAT_INITIAL_MATRIX,&Pt);CHKERRQ(ierr);
      if (viamatmatmatmult) {
        ierr = MatMatMatMult(Pt,A,P,scall,fill,C);CHKERRQ(ierr);
      } else {
        Mat AP;
        ierr = MatMatMult(A,P,MAT_INITIAL_MATRIX,fill,&AP);CHKERRQ(ierr);
        ierr = MatMatMult(Pt,AP,scall,fill,C);CHKERRQ(ierr);
        ierr = MatDestroy(&AP);CHKERRQ(ierr);
      }
      ierr = MatDestroy(&Pt);CHKERRQ(ierr);
    } else {
      ierr = PetscLogEventBegin(MAT_PtAP,A,P,0,0);CHKERRQ(ierr);
      ierr = PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);CHKERRQ(ierr);
      ierr = (*(*C)->ops->ptapnumeric)(A,P,*C);CHKERRQ(ierr);
      ierr = PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);CHKERRQ(ierr);
      ierr = PetscLogEventEnd(MAT_PtAP,A,P,0,0);CHKERRQ(ierr);
    }
    PetscFunctionReturn(0);
  }

  if (fill == PETSC_DEFAULT || fill == PETSC_DECIDE) fill = 2.0;
  if (fill < 1.0) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Expected fill=%g must be >= 1.0",(double)fill);

  fA = A->ops->ptap;
  fP = P->ops->ptap;
  if (fP == fA) {
    if (!fA) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MatPtAP not supported for A of type %s",((PetscObject)A)->type_name);
    ptap = fA;
  } else {
    /* dispatch based on the type of A and P from their PetscObject's PetscFunctionLists. */
    char ptapname[256];
    ierr = PetscStrcpy(ptapname,"MatPtAP_");CHKERRQ(ierr);
    ierr = PetscStrcat(ptapname,((PetscObject)A)->type_name);CHKERRQ(ierr);
    ierr = PetscStrcat(ptapname,"_");CHKERRQ(ierr);
    ierr = PetscStrcat(ptapname,((PetscObject)P)->type_name);CHKERRQ(ierr);
    ierr = PetscStrcat(ptapname,"_C");CHKERRQ(ierr); /* e.g., ptapname = "MatPtAP_seqdense_seqaij_C" */
    ierr = PetscObjectQueryFunction((PetscObject)P,ptapname,&ptap);CHKERRQ(ierr);
    if (!ptap) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"MatPtAP requires A, %s, to be compatible with P, %s",((PetscObject)A)->type_name,((PetscObject)P)->type_name);
  }

  if (viatranspose || viamatmatmatmult) {
    Mat Pt;
    ierr = MatTranspose(P,MAT_INITIAL_MATRIX,&Pt);CHKERRQ(ierr);
    if (viamatmatmatmult) {
      ierr = MatMatMatMult(Pt,A,P,scall,fill,C);CHKERRQ(ierr);
      ierr = PetscInfo(*C,"MatPtAP via MatMatMatMult\n");CHKERRQ(ierr);
    } else {
      Mat AP;
      ierr = MatMatMult(A,P,MAT_INITIAL_MATRIX,fill,&AP);CHKERRQ(ierr);
      ierr = MatMatMult(Pt,AP,scall,fill,C);CHKERRQ(ierr);
      ierr = MatDestroy(&AP);CHKERRQ(ierr);
      ierr = PetscInfo(*C,"MatPtAP via MatTranspose and MatMatMult\n");CHKERRQ(ierr);
    }
    ierr = MatDestroy(&Pt);CHKERRQ(ierr);
  } else {
    ierr = PetscLogEventBegin(MAT_PtAP,A,P,0,0);CHKERRQ(ierr);
    ierr = (*ptap)(A,P,scall,fill,C);CHKERRQ(ierr);
    ierr = PetscLogEventEnd(MAT_PtAP,A,P,0,0);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}


> On Sep 5, 2016, at 4:19 AM, Mark Adams <mfadams at lbl.gov> wrote:
> 
> Barry, 
> 
> THis code looks to me like it should work, but it does not seem to. What code were you looking at?
> 
>   ierr = PetscOptionsGetBool(NULL,NULL,"-matrap",&rap,NULL);CHKERRQ(ierr);
>   if (scall == MAT_INITIAL_MATRIX) {
>     ierr = PetscLogEventBegin(MAT_PtAPSymbolic,A,P,0,0);CHKERRQ(ierr);
>     if (rap) { /* do R=P^T locally, then C=R*A*P */
>       ierr = MatPtAPSymbolic_MPIAIJ_MPIAIJ(A,P,fill,C);CHKERRQ(ierr);
>     } else {       /* do P^T*A*P */
>       ierr = MatPtAPSymbolic_MPIAIJ_MPIAIJ_ptap(A,P,fill,C);CHKERRQ(ierr);
>     }
>     ierr = PetscLogEventEnd(MAT_PtAPSymbolic,A,P,0,0);CHKERRQ(ierr);
>   }
>   ierr = PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);CHKERRQ(ierr);
>   ierr = (*(*C)->ops->ptapnumeric)(A,P,*C);CHKERRQ(ierr);
>   ierr = PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);CHKERRQ(ierr);
> 
> On Sun, Sep 4, 2016 at 9:05 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> 
>    Mark,
> 
>    Look at the source code for MatPtAP. It depends on the approach to how it is computed what logging takes place. I am not sure the rational for this exact logging is, it would seem better to always log it so perhaps that could be changed.
> 
>   Barry
> 
> 
> > On Sep 4, 2016, at 7:10 PM, Mark Adams <mfadams at lbl.gov> wrote:
> >
> > I don't seem to be seeing RAP (MatPtAP, or whatever) profile data. Did this get dumped?
> >
> > I see that there are two timers in GAMG that I can scrap but if RAP is not timed then I will add that to GAMG profile data.
> >
> > Mark
> 
>