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

 18:   MatTransposeGetMat(R,&Rt);
 19:   MatTransposeMatMatMultNumeric_AIJ_AIJ_AIJ_wHYPRE(Rt,A,P,RAP);
 20:   return(0);
 21: }

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

 31:   /* local sizes of matrices will be checked by the calling subroutines */
 32:   MatTransposeGetMat(R,&Rt);
 33:   PetscObjectTypeCompareAny((PetscObject)Rt,&flg,MATSEQAIJ,MATSEQAIJMKL,MATMPIAIJ,NULL);
 34:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)Rt),PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)Rt)->type_name);
 35:   MatTransposeMatMatMultSymbolic_AIJ_AIJ_AIJ_wHYPRE(Rt,A,P,product->fill,RAP);
 36:   RAP->ops->productnumeric = MatProductNumeric_ABC_Transpose_AIJ_AIJ;
 37:   return(0);
 38: }

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

 45:   if (product->type == MATPRODUCT_ABC) {
 46:     C->ops->productsymbolic = MatProductSymbolic_ABC_Transpose_AIJ_AIJ;
 47:   } else SETERRQ1(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"MatProduct type %s is not supported for Transpose, AIJ and AIJ matrices",MatProductTypes[product->type]);
 48:   return(0);
 49: }
 50: #endif

 52: PetscErrorCode MatMatMatMultSymbolic_MPIAIJ_MPIAIJ_MPIAIJ(Mat A,Mat B,Mat C,PetscReal fill,Mat D)
 53: {
 55:   Mat            BC;
 56:   PetscBool      scalable;
 57:   Mat_Product    *product;

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

 78:   D->ops->matmatmultnumeric = MatMatMatMultNumeric_MPIAIJ_MPIAIJ_MPIAIJ;
 79:   return(0);
 80: }

 82: PetscErrorCode MatMatMatMultNumeric_MPIAIJ_MPIAIJ_MPIAIJ(Mat A,Mat B,Mat C,Mat D)
 83: {
 85:   Mat_Product    *product;
 86:   Mat            BC;

 89:   MatCheckProduct(D,4);
 90:   if (!D->product->data) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Product data empty");
 91:   product = D->product;
 92:   BC = product->Dwork;
 93:   if (!BC->ops->matmultnumeric) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Missing numeric operation");
 94:   (*BC->ops->matmultnumeric)(B,C,BC);
 95:   if (!D->ops->matmultnumeric) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Missing numeric operation");
 96:   (*D->ops->matmultnumeric)(A,BC,D);
 97:   return(0);
 98: }

100: /* ----------------------------------------------------- */
101: PetscErrorCode MatDestroy_MPIAIJ_RARt(void *data)
102: {
104:   Mat_RARt       *rart = (Mat_RARt*)data;

107:   MatDestroy(&rart->Rt);
108:   if (rart->destroy) {
109:     (*rart->destroy)(rart->data);
110:   }
111:   PetscFree(rart);
112:   return(0);
113: }

115: PetscErrorCode MatProductNumeric_RARt_MPIAIJ_MPIAIJ(Mat C)
116: {
118:   Mat_RARt       *rart;
119:   Mat            A,R,Rt;

122:   MatCheckProduct(C,1);
123:   if (!C->product->data) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data empty");
124:   rart = (Mat_RARt*)C->product->data;
125:   A    = C->product->A;
126:   R    = C->product->B;
127:   Rt   = rart->Rt;
128:   MatTranspose(R,MAT_REUSE_MATRIX,&Rt);
129:   if (rart->data) C->product->data = rart->data;
130:   (*C->ops->matmatmultnumeric)(R,A,Rt,C);
131:   C->product->data = rart;
132:   return(0);
133: }

135: PetscErrorCode MatProductSymbolic_RARt_MPIAIJ_MPIAIJ(Mat C)
136: {
138:   Mat            A,R,Rt;
139:   Mat_RARt       *rart;

142:   MatCheckProduct(C,1);
143:   if (C->product->data) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data not empty");
144:   A    = C->product->A;
145:   R    = C->product->B;
146:   MatTranspose(R,MAT_INITIAL_MATRIX,&Rt);
147:   /* product->Dwork is used to store A*Rt in MatMatMatMultSymbolic_MPIAIJ_MPIAIJ_MPIAIJ() */
148:   MatMatMatMultSymbolic_MPIAIJ_MPIAIJ_MPIAIJ(R,A,Rt,C->product->fill,C);
149:   C->ops->productnumeric = MatProductNumeric_RARt_MPIAIJ_MPIAIJ;

151:   /* create a supporting struct */
152:   PetscNew(&rart);
153:   rart->Rt      = Rt;
154:   rart->data    = C->product->data;
155:   rart->destroy = C->product->destroy;
156:   C->product->data    = rart;
157:   C->product->destroy = MatDestroy_MPIAIJ_RARt;
158:   return(0);
159: }