1 #define PETSCMAT_DLL 2 3 /* 4 Defines matrix-matrix product routines for pairs of MPIAIJ matrices 5 C = A * B 6 */ 7 #include "src/mat/impls/aij/seq/aij.h" /*I "petscmat.h" I*/ 8 #include "src/mat/utils/freespace.h" 9 #include "src/mat/impls/aij/mpi/mpiaij.h" 10 #include "petscbt.h" 11 12 #undef __FUNCT__ 13 #define __FUNCT__ "MatMatMult_MPIAIJ_MPIAIJ" 14 PetscErrorCode MatMatMult_MPIAIJ_MPIAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill, Mat *C) 15 { 16 PetscErrorCode ierr; 17 18 PetscFunctionBegin; 19 if (scall == MAT_INITIAL_MATRIX){ 20 ierr = MatMatMultSymbolic_MPIAIJ_MPIAIJ(A,B,fill,C);CHKERRQ(ierr);/* numeric product is computed as well */ 21 } else if (scall == MAT_REUSE_MATRIX){ 22 ierr = MatMatMultNumeric_MPIAIJ_MPIAIJ(A,B,*C);CHKERRQ(ierr); 23 } else { 24 SETERRQ1(PETSC_ERR_ARG_WRONG,"Invalid MatReuse %d",(int)scall); 25 } 26 PetscFunctionReturn(0); 27 } 28 29 #undef __FUNCT__ 30 #define __FUNCT__ "PetscObjectContainerDestroy_Mat_MatMatMultMPI" 31 PetscErrorCode PetscObjectContainerDestroy_Mat_MatMatMultMPI(void *ptr) 32 { 33 PetscErrorCode ierr; 34 Mat_MatMatMultMPI *mult=(Mat_MatMatMultMPI*)ptr; 35 36 PetscFunctionBegin; 37 ierr = PetscFree(mult->startsj);CHKERRQ(ierr); 38 ierr = PetscFree(mult->bufa);CHKERRQ(ierr); 39 if (mult->isrowa){ierr = ISDestroy(mult->isrowa);CHKERRQ(ierr);} 40 if (mult->isrowb){ierr = ISDestroy(mult->isrowb);CHKERRQ(ierr);} 41 if (mult->iscolb){ierr = ISDestroy(mult->iscolb);CHKERRQ(ierr);} 42 if (mult->C_seq){ierr = MatDestroy(mult->C_seq);CHKERRQ(ierr);} 43 if (mult->A_loc){ierr = MatDestroy(mult->A_loc);CHKERRQ(ierr); } 44 if (mult->B_seq){ierr = MatDestroy(mult->B_seq);CHKERRQ(ierr);} 45 if (mult->B_loc){ierr = MatDestroy(mult->B_loc);CHKERRQ(ierr);} 46 if (mult->B_oth){ierr = MatDestroy(mult->B_oth);CHKERRQ(ierr);} 47 ierr = PetscFree(mult->abi);CHKERRQ(ierr); 48 ierr = PetscFree(mult->abj);CHKERRQ(ierr); 49 ierr = PetscFree(mult);CHKERRQ(ierr); 50 PetscFunctionReturn(0); 51 } 52 53 EXTERN PetscErrorCode MatDestroy_AIJ(Mat); 54 #undef __FUNCT__ 55 #define __FUNCT__ "MatDestroy_MPIAIJ_MatMatMult" 56 PetscErrorCode MatDestroy_MPIAIJ_MatMatMult(Mat A) 57 { 58 PetscErrorCode ierr; 59 PetscObjectContainer container; 60 Mat_MatMatMultMPI *mult=PETSC_NULL; 61 62 PetscFunctionBegin; 63 ierr = PetscObjectQuery((PetscObject)A,"Mat_MatMatMultMPI",(PetscObject *)&container);CHKERRQ(ierr); 64 if (container) { 65 ierr = PetscObjectContainerGetPointer(container,(void **)&mult);CHKERRQ(ierr); 66 } else { 67 SETERRQ(PETSC_ERR_PLIB,"Container does not exit"); 68 } 69 A->ops->destroy = mult->MatDestroy; 70 ierr = PetscObjectCompose((PetscObject)A,"Mat_MatMatMultMPI",0);CHKERRQ(ierr); 71 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 72 ierr = PetscObjectContainerDestroy(container);CHKERRQ(ierr); 73 PetscFunctionReturn(0); 74 } 75 76 #undef __FUNCT__ 77 #define __FUNCT__ "MatDuplicate_MPIAIJ_MatMatMult" 78 PetscErrorCode MatDuplicate_MPIAIJ_MatMatMult(Mat A, MatDuplicateOption op, Mat *M) { 79 PetscErrorCode ierr; 80 Mat_MatMatMultMPI *mult; 81 PetscObjectContainer container; 82 83 PetscFunctionBegin; 84 ierr = PetscObjectQuery((PetscObject)A,"Mat_MatMatMultMPI",(PetscObject *)&container);CHKERRQ(ierr); 85 if (container) { 86 ierr = PetscObjectContainerGetPointer(container,(void **)&mult);CHKERRQ(ierr); 87 } else { 88 SETERRQ(PETSC_ERR_PLIB,"Container does not exit"); 89 } 90 ierr = (*mult->MatDuplicate)(A,op,M);CHKERRQ(ierr); 91 (*M)->ops->destroy = mult->MatDestroy; /* =MatDestroy_MPIAIJ, *M doesn't duplicate A's container! */ 92 (*M)->ops->duplicate = mult->MatDuplicate; /* =MatDuplicate_ MPIAIJ */ 93 PetscFunctionReturn(0); 94 } 95 96 #undef __FUNCT__ 97 #define __FUNCT__ "MatMatMultSymbolic_MPIAIJ_MPIAIJ" 98 PetscErrorCode MatMatMultSymbolic_MPIAIJ_MPIAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 99 { 100 PetscErrorCode ierr; 101 PetscInt start,end; 102 Mat_MatMatMultMPI *mult; 103 PetscObjectContainer container; 104 105 PetscFunctionBegin; 106 if (A->cmap.rstart != B->rmap.rstart || A->cmap.rend != B->rmap.rend){ 107 SETERRQ4(PETSC_ERR_ARG_SIZ,"Matrix local dimensions are incompatible, (%D, %D) != (%D,%D)",A->cmap.rstart,A->cmap.rend,B->rmap.rstart,B->rmap.rend); 108 } 109 ierr = PetscNew(Mat_MatMatMultMPI,&mult);CHKERRQ(ierr); 110 111 /* create a seq matrix B_seq = submatrix of B by taking rows of B that equal to nonzero col of A */ 112 ierr = MatGetBrowsOfAcols(A,B,MAT_INITIAL_MATRIX,&mult->isrowb,&mult->iscolb,&mult->brstart,&mult->B_seq);CHKERRQ(ierr); 113 114 /* create a seq matrix A_seq = submatrix of A by taking all local rows of A */ 115 start = A->rmap.rstart; end = A->rmap.rend; 116 ierr = ISCreateStride(PETSC_COMM_SELF,end-start,start,1,&mult->isrowa);CHKERRQ(ierr); 117 ierr = MatGetLocalMatCondensed(A,MAT_INITIAL_MATRIX,&mult->isrowa,&mult->isrowb,&mult->A_loc);CHKERRQ(ierr); 118 119 /* compute C_seq = A_seq * B_seq */ 120 ierr = MatMatMult_SeqAIJ_SeqAIJ(mult->A_loc,mult->B_seq,MAT_INITIAL_MATRIX,fill,&mult->C_seq);CHKERRQ(ierr); 121 122 /* create mpi matrix C by concatinating C_seq */ 123 ierr = PetscObjectReference((PetscObject)mult->C_seq);CHKERRQ(ierr); /* prevent C_seq being destroyed by MatMerge() */ 124 ierr = MatMerge(A->comm,mult->C_seq,B->cmap.n,MAT_INITIAL_MATRIX,C);CHKERRQ(ierr); 125 126 /* attach the supporting struct to C for reuse of symbolic C */ 127 ierr = PetscObjectContainerCreate(PETSC_COMM_SELF,&container);CHKERRQ(ierr); 128 ierr = PetscObjectContainerSetPointer(container,mult);CHKERRQ(ierr); 129 ierr = PetscObjectCompose((PetscObject)(*C),"Mat_MatMatMultMPI",(PetscObject)container);CHKERRQ(ierr); 130 ierr = PetscObjectContainerSetUserDestroy(container,PetscObjectContainerDestroy_Mat_MatMatMultMPI);CHKERRQ(ierr); 131 mult->MatDestroy = (*C)->ops->destroy; 132 mult->MatDuplicate = (*C)->ops->duplicate; 133 134 (*C)->ops->destroy = MatDestroy_MPIAIJ_MatMatMult; 135 (*C)->ops->duplicate = MatDuplicate_MPIAIJ_MatMatMult; 136 PetscFunctionReturn(0); 137 } 138 139 /* This routine is called ONLY in the case of reusing previously computed symbolic C */ 140 #undef __FUNCT__ 141 #define __FUNCT__ "MatMatMultNumeric_MPIAIJ_MPIAIJ" 142 PetscErrorCode MatMatMultNumeric_MPIAIJ_MPIAIJ(Mat A,Mat B,Mat C) 143 { 144 PetscErrorCode ierr; 145 Mat *seq; 146 Mat_MatMatMultMPI *mult; 147 PetscObjectContainer container; 148 149 PetscFunctionBegin; 150 ierr = PetscObjectQuery((PetscObject)C,"Mat_MatMatMultMPI",(PetscObject *)&container);CHKERRQ(ierr); 151 if (container) { 152 ierr = PetscObjectContainerGetPointer(container,(void **)&mult);CHKERRQ(ierr); 153 } else { 154 SETERRQ(PETSC_ERR_PLIB,"Container does not exit"); 155 } 156 157 seq = &mult->B_seq; 158 ierr = MatGetSubMatrices(B,1,&mult->isrowb,&mult->iscolb,MAT_REUSE_MATRIX,&seq);CHKERRQ(ierr); 159 mult->B_seq = *seq; 160 161 seq = &mult->A_loc; 162 ierr = MatGetSubMatrices(A,1,&mult->isrowa,&mult->isrowb,MAT_REUSE_MATRIX,&seq);CHKERRQ(ierr); 163 mult->A_loc = *seq; 164 165 ierr = MatMatMult_SeqAIJ_SeqAIJ(mult->A_loc,mult->B_seq,MAT_REUSE_MATRIX,0.0,&mult->C_seq);CHKERRQ(ierr); 166 167 ierr = PetscObjectReference((PetscObject)mult->C_seq);CHKERRQ(ierr); 168 ierr = MatMerge(A->comm,mult->C_seq,B->cmap.n,MAT_REUSE_MATRIX,&C);CHKERRQ(ierr); 169 PetscFunctionReturn(0); 170 } 171