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 MatFreeIntermediateDataStructures_MPIAIJ_BC(Mat ABC)
 34: {
 35:   Mat_MPIAIJ        *a = (Mat_MPIAIJ*)ABC->data;
 36:   Mat_MatMatMatMult *matmatmatmult = a->matmatmatmult;
 37:   PetscErrorCode    ierr;

 40:   if (!matmatmatmult) return(0);

 42:   MatDestroy(&matmatmatmult->BC);
 43:   ABC->ops->destroy = matmatmatmult->destroy;
 44:   PetscFree(a->matmatmatmult);
 45:   return(0);
 46: }

 48: PetscErrorCode MatDestroy_MPIAIJ_MatMatMatMult(Mat A)
 49: {
 50:   PetscErrorCode    ierr;

 53:   (*A->ops->freeintermediatedatastructures)(A);
 54:   (*A->ops->destroy)(A);
 55:   return(0);
 56: }

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

 63:   if (scall == MAT_INITIAL_MATRIX) {
 64:     PetscLogEventBegin(MAT_MatMatMultSymbolic,A,B,C,0);
 65:     MatMatMatMultSymbolic_MPIAIJ_MPIAIJ_MPIAIJ(A,B,C,fill,D);
 66:     PetscLogEventEnd(MAT_MatMatMultSymbolic,A,B,C,0);
 67:   }
 68:   PetscLogEventBegin(MAT_MatMatMultNumeric,A,B,C,0);
 69:   ((*D)->ops->matmatmultnumeric)(A,B,C,*D);
 70:   PetscLogEventEnd(MAT_MatMatMultNumeric,A,B,C,0);
 71:   return(0);
 72: }

 74: PetscErrorCode MatMatMatMultSymbolic_MPIAIJ_MPIAIJ_MPIAIJ(Mat A,Mat B,Mat C,PetscReal fill,Mat *D)
 75: {
 76:   PetscErrorCode    ierr;
 77:   Mat               BC;
 78:   Mat_MatMatMatMult *matmatmatmult;
 79:   Mat_MPIAIJ        *d;
 80:   PetscBool         flg;
 81:   const char        *algTypes[2] = {"scalable","nonscalable"};
 82:   PetscInt          nalg=2,alg=0; /* set default algorithm */

 85:   PetscOptionsBegin(PetscObjectComm((PetscObject)B),((PetscObject)B)->prefix,"MatMatMatMult","Mat");
 86:   PetscOptionsEList("-matmatmatmult_via","Algorithmic approach","MatMatMatMult",algTypes,nalg,algTypes[alg],&alg,&flg);
 87:   PetscOptionsEnd();

 89:   switch (alg) {
 90:   case 0: /* scalable */
 91:     MatMatMultSymbolic_MPIAIJ_MPIAIJ(B,C,fill,&BC);
 92:     MatMatMultSymbolic_MPIAIJ_MPIAIJ(A,BC,fill,D);
 93:     break;
 94:   case 1:
 95:     MatMatMultSymbolic_MPIAIJ_MPIAIJ_nonscalable(B,C,fill,&BC);
 96:     MatMatMultSymbolic_MPIAIJ_MPIAIJ_nonscalable(A,BC,fill,D);
 97:     break;
 98:   default:
 99:     SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Not supported");
100:   }

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

105:   matmatmatmult->BC      = BC;
106:   matmatmatmult->destroy = (*D)->ops->destroy;
107:   d                      = (Mat_MPIAIJ*)(*D)->data;
108:   d->matmatmatmult       = matmatmatmult;

110:   (*D)->ops->matmatmultnumeric = MatMatMatMultNumeric_MPIAIJ_MPIAIJ_MPIAIJ;
111:   (*D)->ops->destroy           = MatDestroy_MPIAIJ_MatMatMatMult;
112:   (*D)->ops->freeintermediatedatastructures = MatFreeIntermediateDataStructures_MPIAIJ_BC;
113:   return(0);
114: }

116: PetscErrorCode MatMatMatMultNumeric_MPIAIJ_MPIAIJ_MPIAIJ(Mat A,Mat B,Mat C,Mat D)
117: {
118:   PetscErrorCode    ierr;
119:   Mat_MPIAIJ        *d             = (Mat_MPIAIJ*)D->data;
120:   Mat_MatMatMatMult *matmatmatmult = d->matmatmatmult;
121:   Mat               BC             = matmatmatmult->BC;

124:   (BC->ops->matmultnumeric)(B,C,BC);
125:   (D->ops->matmultnumeric)(A,BC,D);
126:   return(0);
127: }

129: PetscErrorCode MatRARt_MPIAIJ_MPIAIJ(Mat A,Mat R,MatReuse scall,PetscReal fill,Mat *C)
130: {
132:   Mat            Rt;

135:   MatTranspose(R,MAT_INITIAL_MATRIX,&Rt);
136:   if (scall == MAT_INITIAL_MATRIX) {
137:     MatMatMatMultSymbolic_MPIAIJ_MPIAIJ_MPIAIJ(R,A,Rt,fill,C);
138:   }
139:   MatMatMatMultNumeric_MPIAIJ_MPIAIJ_MPIAIJ(R,A,Rt,*C);
140:   MatDestroy(&Rt);
141:   return(0);
142: }