Actual source code: mpimatmatmult.c

      1: #define PETSCMAT_DLL

      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
     8:  #include src/mat/utils/freespace.h
     9:  #include src/mat/impls/aij/mpi/mpiaij.h
     10:  #include petscbt.h

     14: PetscErrorCode MatMatMult_MPIAIJ_MPIAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill, Mat *C)
     15: {

     19:   if (scall == MAT_INITIAL_MATRIX){
     20:     MatMatMultSymbolic_MPIAIJ_MPIAIJ(A,B,fill,C);/* numeric product is computed as well */
     21:   } else if (scall == MAT_REUSE_MATRIX){
     22:     MatMatMultNumeric_MPIAIJ_MPIAIJ(A,B,*C);
     23:   } else {
     24:     SETERRQ1(PETSC_ERR_ARG_WRONG,"Invalid MatReuse %d",(int)scall);
     25:   }
     26:   return(0);
     27: }

     31: PetscErrorCode PetscObjectContainerDestroy_Mat_MatMatMultMPI(void *ptr)
     32: {
     33:   PetscErrorCode       ierr;
     34:   Mat_MatMatMultMPI    *mult=(Mat_MatMatMultMPI*)ptr;

     37:   if (mult->startsj){PetscFree(mult->startsj);}
     38:   if (mult->bufa){PetscFree(mult->bufa);}
     39:   if (mult->isrowa){ISDestroy(mult->isrowa);}
     40:   if (mult->isrowb){ISDestroy(mult->isrowb);}
     41:   if (mult->iscolb){ISDestroy(mult->iscolb);}
     42:   if (mult->C_seq){MatDestroy(mult->C_seq);}
     43:   if (mult->A_loc){MatDestroy(mult->A_loc); }
     44:   if (mult->B_seq){MatDestroy(mult->B_seq);}
     45:   if (mult->B_loc){MatDestroy(mult->B_loc);}
     46:   if (mult->B_oth){MatDestroy(mult->B_oth);}
     47:   if (mult->abi){PetscFree(mult->abi);}
     48:   if (mult->abj){PetscFree(mult->abj);}
     49:   PetscFree(mult);
     50:   return(0);
     51: }

     53: EXTERN PetscErrorCode MatDestroy_AIJ(Mat);
     56: PetscErrorCode MatDestroy_MPIAIJ_MatMatMult(Mat A)
     57: {
     58:   PetscErrorCode       ierr;
     59:   PetscObjectContainer container;
     60:   Mat_MatMatMultMPI    *mult=PETSC_NULL;

     63:   PetscObjectQuery((PetscObject)A,"Mat_MatMatMultMPI",(PetscObject *)&container);
     64:   if (container) {
     65:     PetscObjectContainerGetPointer(container,(void **)&mult);
     66:   } else {
     67:     SETERRQ(PETSC_ERR_ARG_NULL,"Container does not exit");
     68:   }
     69:   A->ops->destroy = mult->MatDestroy;
     70:   PetscObjectCompose((PetscObject)A,"Mat_MatMatMultMPI",0);
     71:   MatDestroy(A);
     72:   PetscObjectContainerDestroy(container);
     73:   return(0);
     74: }

     78: PetscErrorCode MatMatMultSymbolic_MPIAIJ_MPIAIJ(Mat A,Mat B,PetscReal fill,Mat *C)
     79: {
     80:   Mat_MPIAIJ           *a=(Mat_MPIAIJ*)A->data,*b=(Mat_MPIAIJ*)B->data;
     81:   PetscErrorCode       ierr;
     82:   PetscInt             start,end;
     83:   Mat_MatMatMultMPI    *mult;
     84:   PetscObjectContainer container;
     85:
     87:   if (a->cstart != b->rstart || a->cend != b->rend){
     88:     SETERRQ4(PETSC_ERR_ARG_SIZ,"Matrix local dimensions are incompatible, (%D, %D) != (%D,%D)",a->cstart,a->cend,b->rstart,b->rend);
     89:   }
     90:   PetscNew(Mat_MatMatMultMPI,&mult);

     92:   /* create a seq matrix B_seq = submatrix of B by taking rows of B that equal to nonzero col of A */
     93:   MatGetBrowsOfAcols(A,B,MAT_INITIAL_MATRIX,&mult->isrowb,&mult->iscolb,&mult->brstart,&mult->B_seq);

     95:   /*  create a seq matrix A_seq = submatrix of A by taking all local rows of A */
     96:   start = a->rstart; end = a->rend;
     97:   ISCreateStride(PETSC_COMM_SELF,end-start,start,1,&mult->isrowa);
     98:   MatGetLocalMatCondensed(A,MAT_INITIAL_MATRIX,&mult->isrowa,&mult->isrowb,&mult->A_loc);

    100:   /* compute C_seq = A_seq * B_seq */
    101:   MatMatMult_SeqAIJ_SeqAIJ(mult->A_loc,mult->B_seq,MAT_INITIAL_MATRIX,fill,&mult->C_seq);

    103:   /* create mpi matrix C by concatinating C_seq */
    104:   PetscObjectReference((PetscObject)mult->C_seq); /* prevent C_seq being destroyed by MatMerge() */
    105:   MatMerge(A->comm,mult->C_seq,B->n,MAT_INITIAL_MATRIX,C);
    106:
    107:   /* attach the supporting struct to C for reuse of symbolic C */
    108:   PetscObjectContainerCreate(PETSC_COMM_SELF,&container);
    109:   PetscObjectContainerSetPointer(container,mult);
    110:   PetscObjectCompose((PetscObject)(*C),"Mat_MatMatMultMPI",(PetscObject)container);
    111:   PetscObjectContainerSetUserDestroy(container,PetscObjectContainerDestroy_Mat_MatMatMultMPI);
    112:   mult->MatDestroy = (*C)->ops->destroy;

    114:   (*C)->ops->destroy  = MatDestroy_MPIAIJ_MatMatMult;
    115:   return(0);
    116: }

    118: /* This routine is called ONLY in the case of reusing previously computed symbolic C */
    121: PetscErrorCode MatMatMultNumeric_MPIAIJ_MPIAIJ(Mat A,Mat B,Mat C)
    122: {
    123:   PetscErrorCode       ierr;
    124:   Mat                  *seq;
    125:   Mat_MatMatMultMPI    *mult;
    126:   PetscObjectContainer container;

    129:   PetscObjectQuery((PetscObject)C,"Mat_MatMatMultMPI",(PetscObject *)&container);
    130:   if (container) {
    131:     PetscObjectContainerGetPointer(container,(void **)&mult);
    132:   }

    134:   seq = &mult->B_seq;
    135:   MatGetSubMatrices(B,1,&mult->isrowb,&mult->iscolb,MAT_REUSE_MATRIX,&seq);
    136:   mult->B_seq = *seq;
    137:
    138:   seq = &mult->A_loc;
    139:   MatGetSubMatrices(A,1,&mult->isrowa,&mult->isrowb,MAT_REUSE_MATRIX,&seq);
    140:   mult->A_loc = *seq;

    142:   MatMatMult_SeqAIJ_SeqAIJ(mult->A_loc,mult->B_seq,MAT_REUSE_MATRIX,0.0,&mult->C_seq);

    144:   PetscObjectReference((PetscObject)mult->C_seq);
    145:   MatMerge(A->comm,mult->C_seq,B->n,MAT_REUSE_MATRIX,&C);
    146:   return(0);
    147: }