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(0); 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(0); 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(0); 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,4); 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(0); 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(0); 92 } 93 94 /* ----------------------------------------------------- */ 95 PetscErrorCode MatDestroy_MPIAIJ_RARt(void *data) 96 { 97 Mat_RARt *rart = (Mat_RARt*)data; 98 99 PetscFunctionBegin; 100 PetscCall(MatDestroy(&rart->Rt)); 101 if (rart->destroy) PetscCall((*rart->destroy)(rart->data)); 102 PetscCall(PetscFree(rart)); 103 PetscFunctionReturn(0); 104 } 105 106 PetscErrorCode MatProductNumeric_RARt_MPIAIJ_MPIAIJ(Mat C) 107 { 108 Mat_RARt *rart; 109 Mat A,R,Rt; 110 111 PetscFunctionBegin; 112 MatCheckProduct(C,1); 113 PetscCheck(C->product->data,PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data empty"); 114 rart = (Mat_RARt*)C->product->data; 115 A = C->product->A; 116 R = C->product->B; 117 Rt = rart->Rt; 118 PetscCall(MatTranspose(R,MAT_REUSE_MATRIX,&Rt)); 119 if (rart->data) C->product->data = rart->data; 120 PetscCall((*C->ops->matmatmultnumeric)(R,A,Rt,C)); 121 C->product->data = rart; 122 PetscFunctionReturn(0); 123 } 124 125 PetscErrorCode MatProductSymbolic_RARt_MPIAIJ_MPIAIJ(Mat C) 126 { 127 Mat A,R,Rt; 128 Mat_RARt *rart; 129 130 PetscFunctionBegin; 131 MatCheckProduct(C,1); 132 PetscCheck(!C->product->data,PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data not empty"); 133 A = C->product->A; 134 R = C->product->B; 135 PetscCall(MatTranspose(R,MAT_INITIAL_MATRIX,&Rt)); 136 /* product->Dwork is used to store A*Rt in MatMatMatMultSymbolic_MPIAIJ_MPIAIJ_MPIAIJ() */ 137 PetscCall(MatMatMatMultSymbolic_MPIAIJ_MPIAIJ_MPIAIJ(R,A,Rt,C->product->fill,C)); 138 C->ops->productnumeric = MatProductNumeric_RARt_MPIAIJ_MPIAIJ; 139 140 /* create a supporting struct */ 141 PetscCall(PetscNew(&rart)); 142 rart->Rt = Rt; 143 rart->data = C->product->data; 144 rart->destroy = C->product->destroy; 145 C->product->data = rart; 146 C->product->destroy = MatDestroy_MPIAIJ_RARt; 147 PetscFunctionReturn(0); 148 } 149