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 PetscErrorCode ierr; 14 Mat_Product *product = RAP->product; 15 Mat Rt,R=product->A,A=product->B,P=product->C; 16 17 PetscFunctionBegin; 18 ierr = MatTransposeGetMat(R,&Rt);CHKERRQ(ierr); 19 ierr = MatTransposeMatMatMultNumeric_AIJ_AIJ_AIJ_wHYPRE(Rt,A,P,RAP);CHKERRQ(ierr); 20 PetscFunctionReturn(0); 21 } 22 23 PETSC_INTERN PetscErrorCode MatProductSymbolic_ABC_Transpose_AIJ_AIJ(Mat RAP) 24 { 25 PetscErrorCode ierr; 26 Mat_Product *product = RAP->product; 27 Mat Rt,R=product->A,A=product->B,P=product->C; 28 PetscBool flg; 29 30 PetscFunctionBegin; 31 /* local sizes of matrices will be checked by the calling subroutines */ 32 ierr = MatTransposeGetMat(R,&Rt);CHKERRQ(ierr); 33 ierr = PetscObjectTypeCompareAny((PetscObject)Rt,&flg,MATSEQAIJ,MATSEQAIJMKL,MATMPIAIJ,NULL);CHKERRQ(ierr); 34 if (!flg) SETERRQ1(PetscObjectComm((PetscObject)Rt),PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)Rt)->type_name); 35 ierr = MatTransposeMatMatMultSymbolic_AIJ_AIJ_AIJ_wHYPRE(Rt,A,P,product->fill,RAP);CHKERRQ(ierr); 36 RAP->ops->productnumeric = MatProductNumeric_ABC_Transpose_AIJ_AIJ; 37 PetscFunctionReturn(0); 38 } 39 40 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Transpose_AIJ_AIJ(Mat C) 41 { 42 Mat_Product *product = C->product; 43 44 PetscFunctionBegin; 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 PetscFunctionReturn(0); 49 } 50 #endif 51 52 PetscErrorCode MatMatMatMultSymbolic_MPIAIJ_MPIAIJ_MPIAIJ(Mat A,Mat B,Mat C,PetscReal fill,Mat D) 53 { 54 PetscErrorCode ierr; 55 Mat BC; 56 PetscBool scalable; 57 Mat_Product *product; 58 59 PetscFunctionBegin; 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 ierr = MatProductCreate(B,C,NULL,&BC);CHKERRQ(ierr); 64 ierr = MatProductSetType(BC,MATPRODUCT_AB);CHKERRQ(ierr); 65 ierr = PetscStrcmp(product->alg,"scalable",&scalable);CHKERRQ(ierr); 66 if (scalable) { 67 ierr = MatMatMultSymbolic_MPIAIJ_MPIAIJ(B,C,fill,BC);CHKERRQ(ierr); 68 ierr = MatZeroEntries(BC);CHKERRQ(ierr); /* initialize value entries of BC */ 69 ierr = MatMatMultSymbolic_MPIAIJ_MPIAIJ(A,BC,fill,D);CHKERRQ(ierr); 70 } else { 71 ierr = MatMatMultSymbolic_MPIAIJ_MPIAIJ_nonscalable(B,C,fill,BC);CHKERRQ(ierr); 72 ierr = MatZeroEntries(BC);CHKERRQ(ierr); /* initialize value entries of BC */ 73 ierr = MatMatMultSymbolic_MPIAIJ_MPIAIJ_nonscalable(A,BC,fill,D);CHKERRQ(ierr); 74 } 75 ierr = MatDestroy(&product->Dwork);CHKERRQ(ierr); 76 product->Dwork = BC; 77 78 D->ops->matmatmultnumeric = MatMatMatMultNumeric_MPIAIJ_MPIAIJ_MPIAIJ; 79 PetscFunctionReturn(0); 80 } 81 82 PetscErrorCode MatMatMatMultNumeric_MPIAIJ_MPIAIJ_MPIAIJ(Mat A,Mat B,Mat C,Mat D) 83 { 84 PetscErrorCode ierr; 85 Mat_Product *product; 86 Mat BC; 87 88 PetscFunctionBegin; 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 ierr = (*BC->ops->matmultnumeric)(B,C,BC);CHKERRQ(ierr); 95 if (!D->ops->matmultnumeric) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Missing numeric operation"); 96 ierr = (*D->ops->matmultnumeric)(A,BC,D);CHKERRQ(ierr); 97 PetscFunctionReturn(0); 98 } 99 100 /* ----------------------------------------------------- */ 101 PetscErrorCode MatDestroy_MPIAIJ_RARt(void *data) 102 { 103 PetscErrorCode ierr; 104 Mat_RARt *rart = (Mat_RARt*)data; 105 106 PetscFunctionBegin; 107 ierr = MatDestroy(&rart->Rt);CHKERRQ(ierr); 108 if (rart->destroy) { 109 ierr = (*rart->destroy)(rart->data);CHKERRQ(ierr); 110 } 111 ierr = PetscFree(rart);CHKERRQ(ierr); 112 PetscFunctionReturn(0); 113 } 114 115 PetscErrorCode MatProductNumeric_RARt_MPIAIJ_MPIAIJ(Mat C) 116 { 117 PetscErrorCode ierr; 118 Mat_RARt *rart; 119 Mat A,R,Rt; 120 121 PetscFunctionBegin; 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 ierr = MatTranspose(R,MAT_REUSE_MATRIX,&Rt);CHKERRQ(ierr); 129 if (rart->data) C->product->data = rart->data; 130 ierr = (*C->ops->matmatmultnumeric)(R,A,Rt,C);CHKERRQ(ierr); 131 C->product->data = rart; 132 PetscFunctionReturn(0); 133 } 134 135 PetscErrorCode MatProductSymbolic_RARt_MPIAIJ_MPIAIJ(Mat C) 136 { 137 PetscErrorCode ierr; 138 Mat A,R,Rt; 139 Mat_RARt *rart; 140 141 PetscFunctionBegin; 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 ierr = MatTranspose(R,MAT_INITIAL_MATRIX,&Rt);CHKERRQ(ierr); 147 /* product->Dwork is used to store A*Rt in MatMatMatMultSymbolic_MPIAIJ_MPIAIJ_MPIAIJ() */ 148 ierr = MatMatMatMultSymbolic_MPIAIJ_MPIAIJ_MPIAIJ(R,A,Rt,C->product->fill,C);CHKERRQ(ierr); 149 C->ops->productnumeric = MatProductNumeric_RARt_MPIAIJ_MPIAIJ; 150 151 /* create a supporting struct */ 152 ierr = PetscNew(&rart);CHKERRQ(ierr); 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 PetscFunctionReturn(0); 159 } 160