1: /*
  2:   Defines matrix-matrix-matrix product routines for SeqAIJ matrices
  3:           D = A * B * C
  4: */
  5:  #include <../src/mat/impls/aij/seq/aij.h>

  7: PetscErrorCode MatDestroy_SeqAIJ_MatMatMatMult(Mat A)
  8: {
  9:   Mat_SeqAIJ        *a            = (Mat_SeqAIJ*)A->data;
 10:   Mat_MatMatMatMult *matmatmatmult=a->matmatmatmult;
 11:   PetscErrorCode    ierr;

 14:   MatDestroy(&matmatmatmult->BC);
 15:   matmatmatmult->destroy(A);
 16:   PetscFree(matmatmatmult);
 17:   return(0);
 18: }

 20: PETSC_INTERN PetscErrorCode MatMatMatMult_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,MatReuse scall,PetscReal fill,Mat *D)
 21: {

 25:   if (scall == MAT_INITIAL_MATRIX) {
 26:     PetscLogEventBegin(MAT_MatMatMultSymbolic,A,B,C,0);
 27:     MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(A,B,C,fill,D);
 28:     PetscLogEventEnd(MAT_MatMatMultSymbolic,A,B,C,0);
 29:   }
 30:   PetscLogEventBegin(MAT_MatMatMultNumeric,A,B,C,0);
 31:   ((*D)->ops->matmatmultnumeric)(A,B,C,*D);
 32:   PetscLogEventEnd(MAT_MatMatMultNumeric,A,B,C,0);
 33:   return(0);
 34: }

 36: PetscErrorCode MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,PetscReal fill,Mat *D)
 37: {
 38:   PetscErrorCode    ierr;
 39:   Mat               BC;
 40:   Mat_MatMatMatMult *matmatmatmult;
 41:   Mat_SeqAIJ        *d;
 42:   PetscBool         scalable=PETSC_FALSE;

 45:   PetscObjectOptionsBegin((PetscObject)B);
 46:   PetscOptionsBool("-matmatmatmult_scalable","Use a scalable but slower D=A*B*C","",scalable,&scalable,NULL);
 47:   PetscOptionsEnd();
 48:   if (scalable) {
 49:     MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(B,C,fill,&BC);
 50:     MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(A,BC,fill,D);
 51:   } else {
 52:     MatMatMultSymbolic_SeqAIJ_SeqAIJ(B,C,fill,&BC);
 53:     MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,BC,fill,D);
 54:   }

 56:   /* create struct Mat_MatMatMatMult and attached it to *D */
 57:   PetscNew(&matmatmatmult);

 59:   matmatmatmult->BC      = BC;
 60:   matmatmatmult->destroy = (*D)->ops->destroy;
 61:   d                      = (Mat_SeqAIJ*)(*D)->data;
 62:   d->matmatmatmult       = matmatmatmult;

 64:   (*D)->ops->matmatmultnumeric = MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ;
 65:   (*D)->ops->destroy           = MatDestroy_SeqAIJ_MatMatMatMult;
 66:   return(0);
 67: }

 69: PetscErrorCode MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,Mat D)
 70: {
 71:   PetscErrorCode    ierr;
 72:   Mat_SeqAIJ        *d            =(Mat_SeqAIJ*)D->data;
 73:   Mat_MatMatMatMult *matmatmatmult=d->matmatmatmult;
 74:   Mat               BC            = matmatmatmult->BC;

 77:   (BC->ops->matmultnumeric)(B,C,BC);
 78:   (D->ops->matmultnumeric)(A,BC,D);
 79:   return(0);
 80: }