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