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> /*I "petscmat.h" I*/

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

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

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

 29:   if (scall == MAT_INITIAL_MATRIX) {
 30:     PetscLogEventBegin(MAT_MatMatMultSymbolic,A,B,C,0);
 31:     MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(A,B,C,fill,D);
 32:     PetscLogEventEnd(MAT_MatMatMultSymbolic,A,B,C,0);
 33:   }
 34:   PetscLogEventBegin(MAT_MatMatMultNumeric,A,B,C,0);
 35:   MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(A,B,C,*D);
 36:   PetscLogEventEnd(MAT_MatMatMultNumeric,A,B,C,0);
 37:   return(0);
 38: }

 42: PetscErrorCode MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,PetscReal fill,Mat *D)
 43: {
 44:   PetscErrorCode    ierr;
 45:   Mat               BC;
 46:   Mat_MatMatMatMult *matmatmatmult;
 47:   Mat_SeqAIJ        *d;
 48:   PetscBool         scalable=PETSC_TRUE;

 51:   PetscObjectOptionsBegin((PetscObject)B);
 52:   PetscOptionsBool("-matmatmatmult_scalable","Use a scalable but slower D=A*B*C","",scalable,&scalable,NULL);
 53:   PetscOptionsEnd();
 54:   if (scalable) {
 55:     MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(B,C,fill,&BC);
 56:     MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(A,BC,fill,D);
 57:   } else {
 58:     MatMatMultSymbolic_SeqAIJ_SeqAIJ(B,C,fill,&BC);
 59:     MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,BC,fill,D);
 60:   }

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

 65:   matmatmatmult->BC      = BC;
 66:   matmatmatmult->destroy = (*D)->ops->destroy;
 67:   d                      = (Mat_SeqAIJ*)(*D)->data;
 68:   d->matmatmatmult       = matmatmatmult;

 70:   (*D)->ops->matmatmultnumeric = MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ;
 71:   (*D)->ops->destroy           = MatDestroy_SeqAIJ_MatMatMatMult;
 72:   return(0);
 73: }

 77: PetscErrorCode MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,Mat D)
 78: {
 79:   PetscErrorCode    ierr;
 80:   Mat_SeqAIJ        *d            =(Mat_SeqAIJ*)D->data;
 81:   Mat_MatMatMatMult *matmatmatmult=d->matmatmatmult;
 82:   Mat               BC            = matmatmatmult->BC;

 85:   (BC->ops->matmultnumeric)(B,C,BC);
 86:   (D->ops->matmultnumeric)(A,BC,D);
 87:   return(0);
 88: }