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 PetscErrorCode ierr; 43 Mat_Product *product = C->product; 44 45 PetscFunctionBegin; 46 ierr = MatSetType(C,MATAIJ);CHKERRQ(ierr); 47 if (product->type == MATPRODUCT_ABC) { 48 C->ops->productsymbolic = MatProductSymbolic_ABC_Transpose_AIJ_AIJ; 49 } else SETERRQ1(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"MatProduct type %s is not supported for Transpose, AIJ and AIJ matrices",MatProductTypes[product->type]); 50 PetscFunctionReturn(0); 51 } 52 #endif 53 54 PetscErrorCode MatFreeIntermediateDataStructures_MPIAIJ_BC(Mat ABC) 55 { 56 Mat_MPIAIJ *a = (Mat_MPIAIJ*)ABC->data; 57 Mat_MatMatMatMult *matmatmatmult = a->matmatmatmult; 58 PetscErrorCode ierr; 59 60 PetscFunctionBegin; 61 if (!matmatmatmult) PetscFunctionReturn(0); 62 63 ierr = MatDestroy(&matmatmatmult->BC);CHKERRQ(ierr); 64 ABC->ops->destroy = matmatmatmult->destroy; 65 ierr = PetscFree(a->matmatmatmult);CHKERRQ(ierr); 66 PetscFunctionReturn(0); 67 } 68 69 PetscErrorCode MatDestroy_MPIAIJ_MatMatMatMult(Mat A) 70 { 71 PetscErrorCode ierr; 72 73 PetscFunctionBegin; 74 ierr = (*A->ops->freeintermediatedatastructures)(A);CHKERRQ(ierr); 75 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 76 PetscFunctionReturn(0); 77 } 78 79 PetscErrorCode MatMatMatMultSymbolic_MPIAIJ_MPIAIJ_MPIAIJ(Mat A,Mat B,Mat C,PetscReal fill,Mat D) 80 { 81 PetscErrorCode ierr; 82 Mat BC; 83 PetscBool scalable; 84 Mat_Product *product = D->product; 85 86 PetscFunctionBegin; 87 ierr = MatCreate(PetscObjectComm((PetscObject)A),&BC);CHKERRQ(ierr); 88 if (product) { 89 ierr = PetscStrcmp(product->alg,"scalable",&scalable);CHKERRQ(ierr); 90 } else SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_ARG_NULL,"Call MatProductCreate() first"); 91 92 if (scalable) { 93 ierr = MatMatMultSymbolic_MPIAIJ_MPIAIJ(B,C,fill,BC);CHKERRQ(ierr); 94 ierr = MatZeroEntries(BC);CHKERRQ(ierr); /* initialize value entries of BC */ 95 ierr = MatMatMultSymbolic_MPIAIJ_MPIAIJ(A,BC,fill,D);CHKERRQ(ierr); 96 } else { 97 ierr = MatMatMultSymbolic_MPIAIJ_MPIAIJ_nonscalable(B,C,fill,BC);CHKERRQ(ierr); 98 ierr = MatZeroEntries(BC);CHKERRQ(ierr); /* initialize value entries of BC */ 99 ierr = MatMatMultSymbolic_MPIAIJ_MPIAIJ_nonscalable(A,BC,fill,D);CHKERRQ(ierr); 100 } 101 product->Dwork = BC; 102 103 D->ops->matmatmultnumeric = MatMatMatMultNumeric_MPIAIJ_MPIAIJ_MPIAIJ; 104 D->ops->freeintermediatedatastructures = MatFreeIntermediateDataStructures_MPIAIJ_BC; 105 PetscFunctionReturn(0); 106 } 107 108 PetscErrorCode MatMatMatMultNumeric_MPIAIJ_MPIAIJ_MPIAIJ(Mat A,Mat B,Mat C,Mat D) 109 { 110 PetscErrorCode ierr; 111 Mat_Product *product = D->product; 112 Mat BC = product->Dwork; 113 114 PetscFunctionBegin; 115 ierr = (BC->ops->matmultnumeric)(B,C,BC);CHKERRQ(ierr); 116 ierr = (D->ops->matmultnumeric)(A,BC,D);CHKERRQ(ierr); 117 PetscFunctionReturn(0); 118 } 119 120 /* ----------------------------------------------------- */ 121 PetscErrorCode MatDestroy_MPIAIJ_RARt(Mat C) 122 { 123 PetscErrorCode ierr; 124 Mat_MPIAIJ *c = (Mat_MPIAIJ*)C->data; 125 Mat_RARt *rart = c->rart; 126 127 PetscFunctionBegin; 128 ierr = MatDestroy(&rart->Rt);CHKERRQ(ierr); 129 130 C->ops->destroy = rart->destroy; 131 if (C->ops->destroy) { 132 ierr = (*C->ops->destroy)(C);CHKERRQ(ierr); 133 } 134 ierr = PetscFree(rart);CHKERRQ(ierr); 135 PetscFunctionReturn(0); 136 } 137 138 PetscErrorCode MatProductNumeric_RARt_MPIAIJ_MPIAIJ(Mat C) 139 { 140 PetscErrorCode ierr; 141 Mat_MPIAIJ *c = (Mat_MPIAIJ*)C->data; 142 Mat_RARt *rart = c->rart; 143 Mat_Product *product = C->product; 144 Mat A=product->A,R=product->B,Rt=rart->Rt; 145 146 PetscFunctionBegin; 147 ierr = MatTranspose(R,MAT_REUSE_MATRIX,&Rt);CHKERRQ(ierr); 148 ierr = (C->ops->matmatmultnumeric)(R,A,Rt,C);CHKERRQ(ierr); 149 PetscFunctionReturn(0); 150 } 151 152 PetscErrorCode MatProductSymbolic_RARt_MPIAIJ_MPIAIJ(Mat C) 153 { 154 PetscErrorCode ierr; 155 Mat_Product *product = C->product; 156 Mat A=product->A,R=product->B,Rt; 157 PetscReal fill=product->fill; 158 Mat_RARt *rart; 159 Mat_MPIAIJ *c; 160 161 PetscFunctionBegin; 162 ierr = MatTranspose(R,MAT_INITIAL_MATRIX,&Rt);CHKERRQ(ierr); 163 /* product->Dwork is used to store A*Rt in MatMatMatMultSymbolic_MPIAIJ_MPIAIJ_MPIAIJ() */ 164 ierr = MatMatMatMultSymbolic_MPIAIJ_MPIAIJ_MPIAIJ(R,A,Rt,fill,C);CHKERRQ(ierr); 165 C->ops->productnumeric = MatProductNumeric_RARt_MPIAIJ_MPIAIJ; 166 167 /* create a supporting struct */ 168 ierr = PetscNew(&rart);CHKERRQ(ierr); 169 c = (Mat_MPIAIJ*)C->data; 170 c->rart = rart; 171 rart->Rt = Rt; 172 rart->destroy = C->ops->destroy; 173 C->ops->destroy = MatDestroy_MPIAIJ_RARt; 174 PetscFunctionReturn(0); 175 } 176