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