<div dir="ltr">Thanks,<div>I am actually seeing RAP now ... and the GAMG timers seem to cover GAMG setup pretty well.</div></div><div class="gmail_extra"><br><div class="gmail_quote">On Mon, Sep 5, 2016 at 1:20 PM, Barry Smith <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">#undef __FUNCT__<br>
#define __FUNCT__ "MatPtAP"<br>
/*@<br>
   MatPtAP - Creates the matrix product C = P^T * A * P<br>
<br>
   Neighbor-wise Collective on Mat<br>
<br>
   Input Parameters:<br>
+  A - the matrix<br>
.  P - the projection matrix<br>
.  scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX<br>
-  fill - expected fill as ratio of nnz(C)/(nnz(A) + nnz(P)), use PETSC_DEFAULT if you do not have a good estimate<br>
          if the result is a dense matrix this is irrelevent<br>
<br>
   Output Parameters:<br>
.  C - the product matrix<br>
<br>
   Notes:<br>
   C will be created and must be destroyed by the user with MatDestroy().<br>
<br>
   This routine is currently only implemented for pairs of AIJ matrices and classes<br>
   which inherit from AIJ.<br>
<br>
   Level: intermediate<br>
<br>
.seealso: MatPtAPSymbolic(), MatPtAPNumeric(), MatMatMult(), MatRARt()<br>
@*/<br>
PetscErrorCode MatPtAP(Mat A,Mat P,MatReuse scall,PetscReal fill,Mat *C)<br>
{<br>
  PetscErrorCode ierr;<br>
  PetscErrorCode (*fA)(Mat,Mat,MatReuse,<wbr>PetscReal,Mat*);<br>
  PetscErrorCode (*fP)(Mat,Mat,MatReuse,<wbr>PetscReal,Mat*);<br>
  PetscErrorCode (*ptap)(Mat,Mat,MatReuse,<wbr>PetscReal,Mat*)=NULL;<br>
  PetscBool      viatranspose=PETSC_FALSE,<wbr>viamatmatmatmult=PETSC_FALSE;<br>
<br>
  PetscFunctionBegin;<br>
  ierr = PetscOptionsGetBool(((<wbr>PetscObject)A)->options,((<wbr>PetscObject)A)->prefix,"-<wbr>matptap_viatranspose",&<wbr>viatranspose,NULL);CHKERRQ(<wbr>ierr);<br>
  ierr = PetscOptionsGetBool(((<wbr>PetscObject)A)->options,((<wbr>PetscObject)A)->prefix,"-<wbr>matptap_viamatmatmatmult",&<wbr>viamatmatmatmult,NULL);<wbr>CHKERRQ(ierr);<br>
<br>
  PetscValidHeaderSpecific(A,<wbr>MAT_CLASSID,1);<br>
  PetscValidType(A,1);<br>
  MatCheckPreallocated(A,1);<br>
  if (!A->assembled) SETERRQ(PetscObjectComm((<wbr>PetscObject)A),PETSC_ERR_ARG_<wbr>WRONGSTATE,"Not for unassembled matrix");<br>
  if (A->factortype) SETERRQ(PetscObjectComm((<wbr>PetscObject)A),PETSC_ERR_ARG_<wbr>WRONGSTATE,"Not for factored matrix");<br>
  PetscValidHeaderSpecific(P,<wbr>MAT_CLASSID,2);<br>
  PetscValidType(P,2);<br>
  MatCheckPreallocated(P,2);<br>
  if (!P->assembled) SETERRQ(PetscObjectComm((<wbr>PetscObject)A),PETSC_ERR_ARG_<wbr>WRONGSTATE,"Not for unassembled matrix");<br>
  if (P->factortype) SETERRQ(PetscObjectComm((<wbr>PetscObject)A),PETSC_ERR_ARG_<wbr>WRONGSTATE,"Not for factored matrix");<br>
<br>
  if (A->rmap->N!= A->cmap->N) SETERRQ2(PetscObjectComm((<wbr>PetscObject)A),PETSC_ERR_ARG_<wbr>SIZ,"Matrix A must be square, %D != %D",A->rmap->N,A->cmap->N);<br>
  if (P->rmap->N != A->cmap->N) SETERRQ2(PetscObjectComm((<wbr>PetscObject)A),PETSC_ERR_ARG_<wbr>SIZ,"Matrix dimensions are incompatible, %D != %D",P->rmap->N,A->cmap->N);<br>
  if (fill == PETSC_DEFAULT || fill == PETSC_DECIDE) fill = 2.0;<br>
  if (fill < 1.0) SETERRQ1(PetscObjectComm((<wbr>PetscObject)A),PETSC_ERR_ARG_<wbr>SIZ,"Expected fill=%g must be >= 1.0",(double)fill);<br>
<br>
  if (scall == MAT_REUSE_MATRIX) {<br>
    PetscValidPointer(*C,5);<br>
    PetscValidHeaderSpecific(*C,<wbr>MAT_CLASSID,5);<br>
    if (viatranspose || viamatmatmatmult) {<br>
      Mat Pt;<br>
      ierr = MatTranspose(P,MAT_INITIAL_<wbr>MATRIX,&Pt);CHKERRQ(ierr);<br>
      if (viamatmatmatmult) {<br>
        ierr = MatMatMatMult(Pt,A,P,scall,<wbr>fill,C);CHKERRQ(ierr);<br>
      } else {<br>
        Mat AP;<br>
        ierr = MatMatMult(A,P,MAT_INITIAL_<wbr>MATRIX,fill,&AP);CHKERRQ(ierr)<wbr>;<br>
        ierr = MatMatMult(Pt,AP,scall,fill,C)<wbr>;CHKERRQ(ierr);<br>
        ierr = MatDestroy(&AP);CHKERRQ(ierr);<br>
      }<br>
      ierr = MatDestroy(&Pt);CHKERRQ(ierr);<br>
    } else {<br>
      ierr = PetscLogEventBegin(MAT_PtAP,A,<wbr>P,0,0);CHKERRQ(ierr);<br>
<span class="">      ierr = PetscLogEventBegin(MAT_<wbr>PtAPNumeric,A,P,0,0);CHKERRQ(<wbr>ierr);<br>
      ierr = (*(*C)->ops->ptapnumeric)(A,P,<wbr>*C);CHKERRQ(ierr);<br>
      ierr = PetscLogEventEnd(MAT_<wbr>PtAPNumeric,A,P,0,0);CHKERRQ(<wbr>ierr);<br>
</span>      ierr = PetscLogEventEnd(MAT_PtAP,A,P,<wbr>0,0);CHKERRQ(ierr);<br>
    }<br>
    PetscFunctionReturn(0);<br>
  }<br>
<br>
  if (fill == PETSC_DEFAULT || fill == PETSC_DECIDE) fill = 2.0;<br>
  if (fill < 1.0) SETERRQ1(PetscObjectComm((<wbr>PetscObject)A),PETSC_ERR_ARG_<wbr>SIZ,"Expected fill=%g must be >= 1.0",(double)fill);<br>
<br>
  fA = A->ops->ptap;<br>
  fP = P->ops->ptap;<br>
  if (fP == fA) {<br>
    if (!fA) SETERRQ1(PetscObjectComm((<wbr>PetscObject)A),PETSC_ERR_SUP,"<wbr>MatPtAP not supported for A of type %s",((PetscObject)A)->type_<wbr>name);<br>
    ptap = fA;<br>
  } else {<br>
    /* dispatch based on the type of A and P from their PetscObject's PetscFunctionLists. */<br>
    char ptapname[256];<br>
    ierr = PetscStrcpy(ptapname,"MatPtAP_<wbr>");CHKERRQ(ierr);<br>
    ierr = PetscStrcat(ptapname,((<wbr>PetscObject)A)->type_name);<wbr>CHKERRQ(ierr);<br>
    ierr = PetscStrcat(ptapname,"_");<wbr>CHKERRQ(ierr);<br>
    ierr = PetscStrcat(ptapname,((<wbr>PetscObject)P)->type_name);<wbr>CHKERRQ(ierr);<br>
    ierr = PetscStrcat(ptapname,"_C");<wbr>CHKERRQ(ierr); /* e.g., ptapname = "MatPtAP_seqdense_seqaij_C" */<br>
    ierr = PetscObjectQueryFunction((<wbr>PetscObject)P,ptapname,&ptap);<wbr>CHKERRQ(ierr);<br>
    if (!ptap) SETERRQ2(PetscObjectComm((<wbr>PetscObject)A),PETSC_ERR_ARG_<wbr>INCOMP,"MatPtAP requires A, %s, to be compatible with P, %s",((PetscObject)A)->type_<wbr>name,((PetscObject)P)->type_<wbr>name);<br>
  }<br>
<br>
  if (viatranspose || viamatmatmatmult) {<br>
    Mat Pt;<br>
    ierr = MatTranspose(P,MAT_INITIAL_<wbr>MATRIX,&Pt);CHKERRQ(ierr);<br>
    if (viamatmatmatmult) {<br>
      ierr = MatMatMatMult(Pt,A,P,scall,<wbr>fill,C);CHKERRQ(ierr);<br>
      ierr = PetscInfo(*C,"MatPtAP via MatMatMatMult\n");CHKERRQ(<wbr>ierr);<br>
    } else {<br>
      Mat AP;<br>
      ierr = MatMatMult(A,P,MAT_INITIAL_<wbr>MATRIX,fill,&AP);CHKERRQ(ierr)<wbr>;<br>
      ierr = MatMatMult(Pt,AP,scall,fill,C)<wbr>;CHKERRQ(ierr);<br>
      ierr = MatDestroy(&AP);CHKERRQ(ierr);<br>
      ierr = PetscInfo(*C,"MatPtAP via MatTranspose and MatMatMult\n");CHKERRQ(ierr);<br>
    }<br>
    ierr = MatDestroy(&Pt);CHKERRQ(ierr);<br>
  } else {<br>
    ierr = PetscLogEventBegin(MAT_PtAP,A,<wbr>P,0,0);CHKERRQ(ierr);<br>
    ierr = (*ptap)(A,P,scall,fill,C);<wbr>CHKERRQ(ierr);<br>
    ierr = PetscLogEventEnd(MAT_PtAP,A,P,<wbr>0,0);CHKERRQ(ierr);<br>
  }<br>
  PetscFunctionReturn(0);<br>
<div class="HOEnZb"><div class="h5">}<br>
<br>
<br>
> On Sep 5, 2016, at 4:19 AM, Mark Adams <<a href="mailto:mfadams@lbl.gov">mfadams@lbl.gov</a>> wrote:<br>
><br>
> Barry,<br>
><br>
> THis code looks to me like it should work, but it does not seem to. What code were you looking at?<br>
><br>
>   ierr = PetscOptionsGetBool(NULL,NULL,<wbr>"-matrap",&rap,NULL);CHKERRQ(<wbr>ierr);<br>
>   if (scall == MAT_INITIAL_MATRIX) {<br>
>     ierr = PetscLogEventBegin(MAT_<wbr>PtAPSymbolic,A,P,0,0);CHKERRQ(<wbr>ierr);<br>
>     if (rap) { /* do R=P^T locally, then C=R*A*P */<br>
>       ierr = MatPtAPSymbolic_MPIAIJ_MPIAIJ(<wbr>A,P,fill,C);CHKERRQ(ierr);<br>
>     } else {       /* do P^T*A*P */<br>
>       ierr = MatPtAPSymbolic_MPIAIJ_MPIAIJ_<wbr>ptap(A,P,fill,C);CHKERRQ(ierr)<wbr>;<br>
>     }<br>
>     ierr = PetscLogEventEnd(MAT_<wbr>PtAPSymbolic,A,P,0,0);CHKERRQ(<wbr>ierr);<br>
>   }<br>
>   ierr = PetscLogEventBegin(MAT_<wbr>PtAPNumeric,A,P,0,0);CHKERRQ(<wbr>ierr);<br>
>   ierr = (*(*C)->ops->ptapnumeric)(A,P,<wbr>*C);CHKERRQ(ierr);<br>
>   ierr = PetscLogEventEnd(MAT_<wbr>PtAPNumeric,A,P,0,0);CHKERRQ(<wbr>ierr);<br>
><br>
> On Sun, Sep 4, 2016 at 9:05 PM, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>> wrote:<br>
><br>
>    Mark,<br>
><br>
>    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.<br>
><br>
>   Barry<br>
><br>
><br>
> > On Sep 4, 2016, at 7:10 PM, Mark Adams <<a href="mailto:mfadams@lbl.gov">mfadams@lbl.gov</a>> wrote:<br>
> ><br>
> > I don't seem to be seeing RAP (MatPtAP, or whatever) profile data. Did this get dumped?<br>
> ><br>
> > 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.<br>
> ><br>
> > Mark<br>
><br>
><br>
<br>
</div></div></blockquote></div><br></div>