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> /*I "petscmat.h" I*/ 6 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); 10 11 PETSC_INTERN PetscErrorCode MatProductNumeric_ABC_Transpose_AIJ_AIJ(Mat RAP) { 12 Mat_Product *product = RAP->product; 13 Mat Rt, R = product->A, A = product->B, P = product->C; 14 15 PetscFunctionBegin; 16 PetscCall(MatTransposeGetMat(R, &Rt)); 17 PetscCall(MatTransposeMatMatMultNumeric_AIJ_AIJ_AIJ_wHYPRE(Rt, A, P, RAP)); 18 PetscFunctionReturn(0); 19 } 20 21 PETSC_INTERN PetscErrorCode MatProductSymbolic_ABC_Transpose_AIJ_AIJ(Mat RAP) { 22 Mat_Product *product = RAP->product; 23 Mat Rt, R = product->A, A = product->B, P = product->C; 24 PetscBool flg; 25 26 PetscFunctionBegin; 27 /* local sizes of matrices will be checked by the calling subroutines */ 28 PetscCall(MatTransposeGetMat(R, &Rt)); 29 PetscCall(PetscObjectTypeCompareAny((PetscObject)Rt, &flg, MATSEQAIJ, MATSEQAIJMKL, MATMPIAIJ, NULL)); 30 PetscCheck(flg, PetscObjectComm((PetscObject)Rt), PETSC_ERR_SUP, "Not for matrix type %s", ((PetscObject)Rt)->type_name); 31 PetscCall(MatTransposeMatMatMultSymbolic_AIJ_AIJ_AIJ_wHYPRE(Rt, A, P, product->fill, RAP)); 32 RAP->ops->productnumeric = MatProductNumeric_ABC_Transpose_AIJ_AIJ; 33 PetscFunctionReturn(0); 34 } 35 36 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Transpose_AIJ_AIJ(Mat C) { 37 Mat_Product *product = C->product; 38 39 PetscFunctionBegin; 40 if (product->type == MATPRODUCT_ABC) { 41 C->ops->productsymbolic = MatProductSymbolic_ABC_Transpose_AIJ_AIJ; 42 } else SETERRQ(PetscObjectComm((PetscObject)C), PETSC_ERR_SUP, "MatProduct type %s is not supported for Transpose, AIJ and AIJ matrices", MatProductTypes[product->type]); 43 PetscFunctionReturn(0); 44 } 45 #endif 46 47 PetscErrorCode MatMatMatMultSymbolic_MPIAIJ_MPIAIJ_MPIAIJ(Mat A, Mat B, Mat C, PetscReal fill, Mat D) { 48 Mat BC; 49 PetscBool scalable; 50 Mat_Product *product; 51 52 PetscFunctionBegin; 53 MatCheckProduct(D, 4); 54 PetscCheck(!D->product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data not empty"); 55 product = D->product; 56 PetscCall(MatProductCreate(B, C, NULL, &BC)); 57 PetscCall(MatProductSetType(BC, MATPRODUCT_AB)); 58 PetscCall(PetscStrcmp(product->alg, "scalable", &scalable)); 59 if (scalable) { 60 PetscCall(MatMatMultSymbolic_MPIAIJ_MPIAIJ(B, C, fill, BC)); 61 PetscCall(MatZeroEntries(BC)); /* initialize value entries of BC */ 62 PetscCall(MatMatMultSymbolic_MPIAIJ_MPIAIJ(A, BC, fill, D)); 63 } else { 64 PetscCall(MatMatMultSymbolic_MPIAIJ_MPIAIJ_nonscalable(B, C, fill, BC)); 65 PetscCall(MatZeroEntries(BC)); /* initialize value entries of BC */ 66 PetscCall(MatMatMultSymbolic_MPIAIJ_MPIAIJ_nonscalable(A, BC, fill, D)); 67 } 68 PetscCall(MatDestroy(&product->Dwork)); 69 product->Dwork = BC; 70 71 D->ops->matmatmultnumeric = MatMatMatMultNumeric_MPIAIJ_MPIAIJ_MPIAIJ; 72 PetscFunctionReturn(0); 73 } 74 75 PetscErrorCode MatMatMatMultNumeric_MPIAIJ_MPIAIJ_MPIAIJ(Mat A, Mat B, Mat C, Mat D) { 76 Mat_Product *product; 77 Mat BC; 78 79 PetscFunctionBegin; 80 MatCheckProduct(D, 4); 81 PetscCheck(D->product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data empty"); 82 product = D->product; 83 BC = product->Dwork; 84 PetscCall((*BC->ops->matmultnumeric)(B, C, BC)); 85 PetscCall((*D->ops->matmultnumeric)(A, BC, D)); 86 PetscFunctionReturn(0); 87 } 88 89 /* ----------------------------------------------------- */ 90 PetscErrorCode MatDestroy_MPIAIJ_RARt(void *data) { 91 Mat_RARt *rart = (Mat_RARt *)data; 92 93 PetscFunctionBegin; 94 PetscCall(MatDestroy(&rart->Rt)); 95 if (rart->destroy) PetscCall((*rart->destroy)(rart->data)); 96 PetscCall(PetscFree(rart)); 97 PetscFunctionReturn(0); 98 } 99 100 PetscErrorCode MatProductNumeric_RARt_MPIAIJ_MPIAIJ(Mat C) { 101 Mat_RARt *rart; 102 Mat A, R, Rt; 103 104 PetscFunctionBegin; 105 MatCheckProduct(C, 1); 106 PetscCheck(C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data empty"); 107 rart = (Mat_RARt *)C->product->data; 108 A = C->product->A; 109 R = C->product->B; 110 Rt = rart->Rt; 111 PetscCall(MatTranspose(R, MAT_REUSE_MATRIX, &Rt)); 112 if (rart->data) C->product->data = rart->data; 113 PetscCall((*C->ops->matmatmultnumeric)(R, A, Rt, C)); 114 C->product->data = rart; 115 PetscFunctionReturn(0); 116 } 117 118 PetscErrorCode MatProductSymbolic_RARt_MPIAIJ_MPIAIJ(Mat C) { 119 Mat A, R, Rt; 120 Mat_RARt *rart; 121 122 PetscFunctionBegin; 123 MatCheckProduct(C, 1); 124 PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty"); 125 A = C->product->A; 126 R = C->product->B; 127 PetscCall(MatTranspose(R, MAT_INITIAL_MATRIX, &Rt)); 128 /* product->Dwork is used to store A*Rt in MatMatMatMultSymbolic_MPIAIJ_MPIAIJ_MPIAIJ() */ 129 PetscCall(MatMatMatMultSymbolic_MPIAIJ_MPIAIJ_MPIAIJ(R, A, Rt, C->product->fill, C)); 130 C->ops->productnumeric = MatProductNumeric_RARt_MPIAIJ_MPIAIJ; 131 132 /* create a supporting struct */ 133 PetscCall(PetscNew(&rart)); 134 rart->Rt = Rt; 135 rart->data = C->product->data; 136 rart->destroy = C->product->destroy; 137 C->product->data = rart; 138 C->product->destroy = MatDestroy_MPIAIJ_RARt; 139 PetscFunctionReturn(0); 140 } 141