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