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 { 13 Mat_Product *product = RAP->product; 14 Mat Rt, R = product->A, A = product->B, P = product->C; 15 16 PetscFunctionBegin; 17 PetscCall(MatTransposeGetMat(R, &Rt)); 18 PetscCall(MatTransposeMatMatMultNumeric_AIJ_AIJ_AIJ_wHYPRE(Rt, A, P, RAP)); 19 PetscFunctionReturn(PETSC_SUCCESS); 20 } 21 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; 27 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 } 37 38 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Transpose_AIJ_AIJ(Mat C) 39 { 40 Mat_Product *product = C->product; 41 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 49 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; 55 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; 74 75 D->ops->matmatmultnumeric = MatMatMatMultNumeric_MPIAIJ_MPIAIJ_MPIAIJ; 76 PetscFunctionReturn(PETSC_SUCCESS); 77 } 78 79 PetscErrorCode MatMatMatMultNumeric_MPIAIJ_MPIAIJ_MPIAIJ(Mat A, Mat B, Mat C, Mat D) 80 { 81 Mat_Product *product; 82 Mat BC; 83 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 } 93 94 PetscErrorCode MatDestroy_MPIAIJ_RARt(void *data) 95 { 96 Mat_RARt *rart = (Mat_RARt *)data; 97 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 } 104 105 PetscErrorCode MatProductNumeric_RARt_MPIAIJ_MPIAIJ(Mat C) 106 { 107 Mat_RARt *rart; 108 Mat A, R, Rt; 109 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 } 123 124 PetscErrorCode MatProductSymbolic_RARt_MPIAIJ_MPIAIJ(Mat C) 125 { 126 Mat A, R, Rt; 127 Mat_RARt *rart; 128 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; 138 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 } 148