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 MatProductNumeric_ABC_Transpose_AIJ_AIJ(Mat RAP)
 12: {
 13:   Mat_Product *product = RAP->product;
 14:   Mat          Rt, R = product->A, A = product->B, P = product->C;

 16:   PetscFunctionBegin;
 17:   PetscCall(MatTransposeGetMat(R, &Rt));
 18:   PetscCall(MatTransposeMatMatMultNumeric_AIJ_AIJ_AIJ_wHYPRE(Rt, A, P, RAP));
 19:   PetscFunctionReturn(PETSC_SUCCESS);
 20: }

 22: PETSC_INTERN PetscErrorCode MatProductSymbolic_ABC_Transpose_AIJ_AIJ(Mat RAP)
 23: {
 24:   Mat_Product *product = RAP->product;
 25:   Mat          Rt, R = product->A, A = product->B, P = product->C;
 26:   PetscBool    flg;

 28:   PetscFunctionBegin;
 29:   /* local sizes of matrices will be checked by the calling subroutines */
 30:   PetscCall(MatTransposeGetMat(R, &Rt));
 31:   PetscCall(PetscObjectTypeCompareAny((PetscObject)Rt, &flg, MATSEQAIJ, MATSEQAIJMKL, MATMPIAIJ, NULL));
 32:   PetscCheck(flg, PetscObjectComm((PetscObject)Rt), PETSC_ERR_SUP, "Not for matrix type %s", ((PetscObject)Rt)->type_name);
 33:   PetscCall(MatTransposeMatMatMultSymbolic_AIJ_AIJ_AIJ_wHYPRE(Rt, A, P, product->fill, RAP));
 34:   RAP->ops->productnumeric = MatProductNumeric_ABC_Transpose_AIJ_AIJ;
 35:   PetscFunctionReturn(PETSC_SUCCESS);
 36: }

 38: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Transpose_AIJ_AIJ(Mat C)
 39: {
 40:   Mat_Product *product = C->product;

 42:   PetscFunctionBegin;
 43:   if (product->type == MATPRODUCT_ABC) {
 44:     C->ops->productsymbolic = MatProductSymbolic_ABC_Transpose_AIJ_AIJ;
 45:   } else SETERRQ(PetscObjectComm((PetscObject)C), PETSC_ERR_SUP, "MatProduct type %s is not supported for Transpose, AIJ and AIJ matrices", MatProductTypes[product->type]);
 46:   PetscFunctionReturn(PETSC_SUCCESS);
 47: }
 48: #endif

 50: PetscErrorCode MatMatMatMultSymbolic_MPIAIJ_MPIAIJ_MPIAIJ(Mat A, Mat B, Mat C, PetscReal fill, Mat D)
 51: {
 52:   Mat          BC;
 53:   PetscBool    scalable;
 54:   Mat_Product *product;

 56:   PetscFunctionBegin;
 57:   MatCheckProduct(D, 5);
 58:   PetscCheck(!D->product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data not empty");
 59:   product = D->product;
 60:   PetscCall(MatProductCreate(B, C, NULL, &BC));
 61:   PetscCall(MatProductSetType(BC, MATPRODUCT_AB));
 62:   PetscCall(PetscStrcmp(product->alg, "scalable", &scalable));
 63:   if (scalable) {
 64:     PetscCall(MatMatMultSymbolic_MPIAIJ_MPIAIJ(B, C, fill, BC));
 65:     PetscCall(MatZeroEntries(BC)); /* initialize value entries of BC */
 66:     PetscCall(MatMatMultSymbolic_MPIAIJ_MPIAIJ(A, BC, fill, D));
 67:   } else {
 68:     PetscCall(MatMatMultSymbolic_MPIAIJ_MPIAIJ_nonscalable(B, C, fill, BC));
 69:     PetscCall(MatZeroEntries(BC)); /* initialize value entries of BC */
 70:     PetscCall(MatMatMultSymbolic_MPIAIJ_MPIAIJ_nonscalable(A, BC, fill, D));
 71:   }
 72:   PetscCall(MatDestroy(&product->Dwork));
 73:   product->Dwork = BC;

 75:   D->ops->matmatmultnumeric = MatMatMatMultNumeric_MPIAIJ_MPIAIJ_MPIAIJ;
 76:   PetscFunctionReturn(PETSC_SUCCESS);
 77: }

 79: PetscErrorCode MatMatMatMultNumeric_MPIAIJ_MPIAIJ_MPIAIJ(Mat A, Mat B, Mat C, Mat D)
 80: {
 81:   Mat_Product *product;
 82:   Mat          BC;

 84:   PetscFunctionBegin;
 85:   MatCheckProduct(D, 4);
 86:   PetscCheck(D->product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data empty");
 87:   product = D->product;
 88:   BC      = product->Dwork;
 89:   PetscCall((*BC->ops->matmultnumeric)(B, C, BC));
 90:   PetscCall((*D->ops->matmultnumeric)(A, BC, D));
 91:   PetscFunctionReturn(PETSC_SUCCESS);
 92: }

 94: static PetscErrorCode MatDestroy_MPIAIJ_RARt(void *data)
 95: {
 96:   Mat_RARt *rart = (Mat_RARt *)data;

 98:   PetscFunctionBegin;
 99:   PetscCall(MatDestroy(&rart->Rt));
100:   if (rart->destroy) PetscCall((*rart->destroy)(rart->data));
101:   PetscCall(PetscFree(rart));
102:   PetscFunctionReturn(PETSC_SUCCESS);
103: }

105: PetscErrorCode MatProductNumeric_RARt_MPIAIJ_MPIAIJ(Mat C)
106: {
107:   Mat_RARt *rart;
108:   Mat       A, R, Rt;

110:   PetscFunctionBegin;
111:   MatCheckProduct(C, 1);
112:   PetscCheck(C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data empty");
113:   rart = (Mat_RARt *)C->product->data;
114:   A    = C->product->A;
115:   R    = C->product->B;
116:   Rt   = rart->Rt;
117:   PetscCall(MatTranspose(R, MAT_REUSE_MATRIX, &Rt));
118:   if (rart->data) C->product->data = rart->data;
119:   PetscCall((*C->ops->matmatmultnumeric)(R, A, Rt, C));
120:   C->product->data = rart;
121:   PetscFunctionReturn(PETSC_SUCCESS);
122: }

124: PetscErrorCode MatProductSymbolic_RARt_MPIAIJ_MPIAIJ(Mat C)
125: {
126:   Mat       A, R, Rt;
127:   Mat_RARt *rart;

129:   PetscFunctionBegin;
130:   MatCheckProduct(C, 1);
131:   PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
132:   A = C->product->A;
133:   R = C->product->B;
134:   PetscCall(MatTranspose(R, MAT_INITIAL_MATRIX, &Rt));
135:   /* product->Dwork is used to store A*Rt in MatMatMatMultSymbolic_MPIAIJ_MPIAIJ_MPIAIJ() */
136:   PetscCall(MatMatMatMultSymbolic_MPIAIJ_MPIAIJ_MPIAIJ(R, A, Rt, C->product->fill, C));
137:   C->ops->productnumeric = MatProductNumeric_RARt_MPIAIJ_MPIAIJ;

139:   /* create a supporting struct */
140:   PetscCall(PetscNew(&rart));
141:   rart->Rt            = Rt;
142:   rart->data          = C->product->data;
143:   rart->destroy       = C->product->destroy;
144:   C->product->data    = rart;
145:   C->product->destroy = MatDestroy_MPIAIJ_RARt;
146:   PetscFunctionReturn(PETSC_SUCCESS);
147: }