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:     MatZeroEntries(BC); /* initialize value entries of BC */
 93:     MatMatMultSymbolic_MPIAIJ_MPIAIJ(A,BC,fill,D);
 94:     break;
 95:   case 1:
 96:     MatMatMultSymbolic_MPIAIJ_MPIAIJ_nonscalable(B,C,fill,&BC);
 97:     MatZeroEntries(BC); /* initialize value entries of BC */
 98:     MatMatMultSymbolic_MPIAIJ_MPIAIJ_nonscalable(A,BC,fill,D);
 99:     break;
100:   default:
101:     SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Not supported");
102:   }

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

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

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

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

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

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

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