1: /*
  2:   Defines matrix-matrix-matrix product routines for MPIAIJ matrices
  3:           D = A * B * C
  4: */
  5:  #include <../src/mat/impls/aij/mpi/mpiaij.h>

  7: #if defined(PETSC_HAVE_HYPRE)
  8: PETSC_INTERN PetscErrorCode MatTransposeMatMatMultSymbolic_AIJ_AIJ_AIJ_wHYPRE(Mat,Mat,Mat,PetscReal,Mat*);
  9: PETSC_INTERN PetscErrorCode MatTransposeMatMatMultNumeric_AIJ_AIJ_AIJ_wHYPRE(Mat,Mat,Mat,Mat);

 11: PETSC_INTERN PetscErrorCode MatMatMatMult_Transpose_AIJ_AIJ(Mat R,Mat A,Mat P,MatReuse scall, PetscReal fill, Mat *RAP)
 12: {
 13:   Mat            Rt;
 14:   PetscBool      flg;

 18:   MatTransposeGetMat(R,&Rt);
 19:   PetscObjectTypeCompareAny((PetscObject)Rt,&flg,MATSEQAIJ,MATMPIAIJ,NULL);
 20:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)Rt),PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)Rt)->type_name);
 21:   if (scall == MAT_INITIAL_MATRIX) {
 22:     PetscLogEventBegin(MAT_MatMatMultSymbolic,R,A,P,0);
 23:     MatTransposeMatMatMultSymbolic_AIJ_AIJ_AIJ_wHYPRE(Rt,A,P,fill,RAP);
 24:     PetscLogEventEnd(MAT_MatMatMultSymbolic,R,A,P,0);
 25:   }
 26:   PetscLogEventBegin(MAT_MatMatMultNumeric,R,A,P,0);
 27:   MatTransposeMatMatMultNumeric_AIJ_AIJ_AIJ_wHYPRE(Rt,A,P,*RAP);
 28:   PetscLogEventEnd(MAT_MatMatMultNumeric,R,A,P,0);
 29:   return(0);
 30: }
 31: #endif

 33: PetscErrorCode MatDestroy_MPIAIJ_MatMatMatMult(Mat A)
 34: {
 35:   Mat_MPIAIJ        *a            = (Mat_MPIAIJ*)A->data;
 36:   Mat_MatMatMatMult *matmatmatmult=a->matmatmatmult;
 37:   PetscErrorCode    ierr;

 40:   MatDestroy(&matmatmatmult->BC);
 41:   matmatmatmult->destroy(A);
 42:   PetscFree(matmatmatmult);
 43:   return(0);
 44: }

 46: PETSC_INTERN PetscErrorCode MatMatMatMult_MPIAIJ_MPIAIJ_MPIAIJ(Mat A,Mat B,Mat C,MatReuse scall,PetscReal fill,Mat *D)
 47: {

 51:   if (scall == MAT_INITIAL_MATRIX) {
 52:     PetscLogEventBegin(MAT_MatMatMultSymbolic,A,B,C,0);
 53:     MatMatMatMultSymbolic_MPIAIJ_MPIAIJ_MPIAIJ(A,B,C,fill,D);
 54:     PetscLogEventEnd(MAT_MatMatMultSymbolic,A,B,C,0);
 55:   }
 56:   PetscLogEventBegin(MAT_MatMatMultNumeric,A,B,C,0);
 57:   ((*D)->ops->matmatmultnumeric)(A,B,C,*D);
 58:   PetscLogEventEnd(MAT_MatMatMultNumeric,A,B,C,0);
 59:   return(0);
 60: }

 62: PetscErrorCode MatMatMatMultSymbolic_MPIAIJ_MPIAIJ_MPIAIJ(Mat A,Mat B,Mat C,PetscReal fill,Mat *D)
 63: {
 64:   PetscErrorCode    ierr;
 65:   Mat               BC;
 66:   Mat_MatMatMatMult *matmatmatmult;
 67:   Mat_MPIAIJ        *d;
 68:   PetscBool         scalable=PETSC_TRUE;

 71:   PetscObjectOptionsBegin((PetscObject)B);
 72:   PetscOptionsBool("-matmatmatmult_scalable","Use a scalable but slower D=A*B*C","",scalable,&scalable,NULL);
 73:   PetscOptionsEnd();
 74:   if (scalable) {
 75:     MatMatMultSymbolic_MPIAIJ_MPIAIJ(B,C,fill,&BC);
 76:     MatMatMultSymbolic_MPIAIJ_MPIAIJ(A,BC,fill,D);
 77:   } else {
 78:     MatMatMultSymbolic_MPIAIJ_MPIAIJ_nonscalable(B,C,fill,&BC);
 79:     MatMatMultSymbolic_MPIAIJ_MPIAIJ_nonscalable(A,BC,fill,D);
 80:   }

 82:   /* create struct Mat_MatMatMatMult and attached it to *D */
 83:   PetscNew(&matmatmatmult);

 85:   matmatmatmult->BC      = BC;
 86:   matmatmatmult->destroy = (*D)->ops->destroy;
 87:   d                      = (Mat_MPIAIJ*)(*D)->data;
 88:   d->matmatmatmult       = matmatmatmult;

 90:   (*D)->ops->matmatmultnumeric = MatMatMatMultNumeric_MPIAIJ_MPIAIJ_MPIAIJ;
 91:   (*D)->ops->destroy           = MatDestroy_MPIAIJ_MatMatMatMult;
 92:   return(0);
 93: }

 95: PetscErrorCode MatMatMatMultNumeric_MPIAIJ_MPIAIJ_MPIAIJ(Mat A,Mat B,Mat C,Mat D)
 96: {
 97:   PetscErrorCode    ierr;
 98:   Mat_MPIAIJ        *d             = (Mat_MPIAIJ*)D->data;
 99:   Mat_MatMatMatMult *matmatmatmult = d->matmatmatmult;
100:   Mat               BC             = matmatmatmult->BC;

103:   (BC->ops->matmultnumeric)(B,C,BC);
104:   (D->ops->matmultnumeric)(A,BC,D);
105:   return(0);
106: }

108: PetscErrorCode MatRARt_MPIAIJ_MPIAIJ(Mat A,Mat R,MatReuse scall,PetscReal fill,Mat *C)
109: {
111:   Mat            Rt;

114:   MatTranspose(R,MAT_INITIAL_MATRIX,&Rt);
115:   if (scall == MAT_INITIAL_MATRIX) {
116:     MatMatMatMultSymbolic_MPIAIJ_MPIAIJ_MPIAIJ(R,A,Rt,fill,C);
117:   }
118:   MatMatMatMultNumeric_MPIAIJ_MPIAIJ_MPIAIJ(R,A,Rt,*C);
119:   MatDestroy(&Rt);
120:   return(0);
121: }