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 MatMatMatMult_Transpose_AIJ_AIJ(Mat R,Mat A,Mat P,MatReuse scall, PetscReal fill, Mat *RAP) 12 { 13 Mat Rt; 14 PetscBool flg; 15 PetscErrorCode ierr; 16 17 PetscFunctionBegin; 18 ierr = MatTransposeGetMat(R,&Rt);CHKERRQ(ierr); 19 ierr = PetscObjectTypeCompareAny((PetscObject)Rt,&flg,MATSEQAIJ,MATMPIAIJ,NULL);CHKERRQ(ierr); 20 if (!flg) SETERRQ1(PetscObjectComm((PetscObject)Rt),PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)Rt)->type_name); 21 if (scall == MAT_INITIAL_MATRIX) { 22 ierr = PetscLogEventBegin(MAT_MatMatMultSymbolic,R,A,P,0);CHKERRQ(ierr); 23 ierr = MatTransposeMatMatMultSymbolic_AIJ_AIJ_AIJ_wHYPRE(Rt,A,P,fill,RAP);CHKERRQ(ierr); 24 ierr = PetscLogEventEnd(MAT_MatMatMultSymbolic,R,A,P,0);CHKERRQ(ierr); 25 } 26 ierr = PetscLogEventBegin(MAT_MatMatMultNumeric,R,A,P,0);CHKERRQ(ierr); 27 ierr = MatTransposeMatMatMultNumeric_AIJ_AIJ_AIJ_wHYPRE(Rt,A,P,*RAP);CHKERRQ(ierr); 28 ierr = PetscLogEventEnd(MAT_MatMatMultNumeric,R,A,P,0);CHKERRQ(ierr); 29 PetscFunctionReturn(0); 30 } 31 #endif 32 33 PetscErrorCode MatFreeIntermediateDataStructures_MPIAIJ_BC(Mat ABC) 34 { 35 Mat_MPIAIJ *a = (Mat_MPIAIJ*)ABC->data; 36 Mat_MatMatMatMult *matmatmatmult = a->matmatmatmult; 37 PetscErrorCode ierr; 38 39 PetscFunctionBegin; 40 if (!matmatmatmult) PetscFunctionReturn(0); 41 42 ierr = MatDestroy(&matmatmatmult->BC);CHKERRQ(ierr); 43 ABC->ops->destroy = matmatmatmult->destroy; 44 ierr = PetscFree(a->matmatmatmult);CHKERRQ(ierr); 45 PetscFunctionReturn(0); 46 } 47 48 PetscErrorCode MatDestroy_MPIAIJ_MatMatMatMult(Mat A) 49 { 50 PetscErrorCode ierr; 51 52 PetscFunctionBegin; 53 ierr = (*A->ops->freeintermediatedatastructures)(A);CHKERRQ(ierr); 54 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 55 PetscFunctionReturn(0); 56 } 57 58 PETSC_INTERN PetscErrorCode MatMatMatMult_MPIAIJ_MPIAIJ_MPIAIJ(Mat A,Mat B,Mat C,MatReuse scall,PetscReal fill,Mat *D) 59 { 60 PetscErrorCode ierr; 61 62 PetscFunctionBegin; 63 if (scall == MAT_INITIAL_MATRIX) { 64 ierr = PetscLogEventBegin(MAT_MatMatMultSymbolic,A,B,C,0);CHKERRQ(ierr); 65 ierr = MatMatMatMultSymbolic_MPIAIJ_MPIAIJ_MPIAIJ(A,B,C,fill,D);CHKERRQ(ierr); 66 ierr = PetscLogEventEnd(MAT_MatMatMultSymbolic,A,B,C,0);CHKERRQ(ierr); 67 } 68 ierr = PetscLogEventBegin(MAT_MatMatMultNumeric,A,B,C,0);CHKERRQ(ierr); 69 ierr = ((*D)->ops->matmatmultnumeric)(A,B,C,*D);CHKERRQ(ierr); 70 ierr = PetscLogEventEnd(MAT_MatMatMultNumeric,A,B,C,0);CHKERRQ(ierr); 71 PetscFunctionReturn(0); 72 } 73 74 PetscErrorCode MatMatMatMultSymbolic_MPIAIJ_MPIAIJ_MPIAIJ(Mat A,Mat B,Mat C,PetscReal fill,Mat *D) 75 { 76 PetscErrorCode ierr; 77 Mat BC; 78 Mat_MatMatMatMult *matmatmatmult; 79 Mat_MPIAIJ *d; 80 PetscBool flg; 81 const char *algTypes[2] = {"scalable","nonscalable"}; 82 PetscInt nalg=2,alg=0; /* set default algorithm */ 83 84 PetscFunctionBegin; 85 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)B),((PetscObject)B)->prefix,"MatMatMatMult","Mat");CHKERRQ(ierr); 86 ierr = PetscOptionsEList("-matmatmatmult_via","Algorithmic approach","MatMatMatMult",algTypes,nalg,algTypes[alg],&alg,&flg);CHKERRQ(ierr); 87 ierr = PetscOptionsEnd();CHKERRQ(ierr); 88 89 switch (alg) { 90 case 0: /* scalable */ 91 ierr = MatMatMultSymbolic_MPIAIJ_MPIAIJ(B,C,fill,&BC);CHKERRQ(ierr); 92 ierr = MatMatMultSymbolic_MPIAIJ_MPIAIJ(A,BC,fill,D);CHKERRQ(ierr); 93 break; 94 case 1: 95 ierr = MatMatMultSymbolic_MPIAIJ_MPIAIJ_nonscalable(B,C,fill,&BC);CHKERRQ(ierr); 96 ierr = MatMatMultSymbolic_MPIAIJ_MPIAIJ_nonscalable(A,BC,fill,D);CHKERRQ(ierr); 97 break; 98 default: 99 SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Not supported"); 100 } 101 102 /* create struct Mat_MatMatMatMult and attached it to *D */ 103 ierr = PetscNew(&matmatmatmult);CHKERRQ(ierr); 104 105 matmatmatmult->BC = BC; 106 matmatmatmult->destroy = (*D)->ops->destroy; 107 d = (Mat_MPIAIJ*)(*D)->data; 108 d->matmatmatmult = matmatmatmult; 109 110 (*D)->ops->matmatmultnumeric = MatMatMatMultNumeric_MPIAIJ_MPIAIJ_MPIAIJ; 111 (*D)->ops->destroy = MatDestroy_MPIAIJ_MatMatMatMult; 112 (*D)->ops->freeintermediatedatastructures = MatFreeIntermediateDataStructures_MPIAIJ_BC; 113 PetscFunctionReturn(0); 114 } 115 116 PetscErrorCode MatMatMatMultNumeric_MPIAIJ_MPIAIJ_MPIAIJ(Mat A,Mat B,Mat C,Mat D) 117 { 118 PetscErrorCode ierr; 119 Mat_MPIAIJ *d = (Mat_MPIAIJ*)D->data; 120 Mat_MatMatMatMult *matmatmatmult = d->matmatmatmult; 121 Mat BC = matmatmatmult->BC; 122 123 PetscFunctionBegin; 124 ierr = (BC->ops->matmultnumeric)(B,C,BC);CHKERRQ(ierr); 125 ierr = (D->ops->matmultnumeric)(A,BC,D);CHKERRQ(ierr); 126 PetscFunctionReturn(0); 127 } 128 129 PetscErrorCode MatRARt_MPIAIJ_MPIAIJ(Mat A,Mat R,MatReuse scall,PetscReal fill,Mat *C) 130 { 131 PetscErrorCode ierr; 132 Mat Rt; 133 134 PetscFunctionBegin; 135 ierr = MatTranspose(R,MAT_INITIAL_MATRIX,&Rt);CHKERRQ(ierr); 136 if (scall == MAT_INITIAL_MATRIX) { 137 ierr = MatMatMatMultSymbolic_MPIAIJ_MPIAIJ_MPIAIJ(R,A,Rt,fill,C);CHKERRQ(ierr); 138 } 139 ierr = MatMatMatMultNumeric_MPIAIJ_MPIAIJ_MPIAIJ(R,A,Rt,*C);CHKERRQ(ierr); 140 ierr = MatDestroy(&Rt);CHKERRQ(ierr); 141 PetscFunctionReturn(0); 142 } 143