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: PetscErrorCode MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,PetscReal fill,Mat D)
 21: {
 22:   PetscErrorCode    ierr;
 23:   Mat               BC;
 24:   Mat_MatMatMatMult *matmatmatmult;
 25:   Mat_SeqAIJ        *d;
 26:   Mat_Product       *product = D->product;
 27:   MatProductAlgorithm alg=product->alg;

 30:   if (!product) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Data struc Mat_Product is not created, call MatProductCreate() first");
 31:   MatCreate(PETSC_COMM_SELF,&BC);
 32:   MatMatMultSymbolic_SeqAIJ_SeqAIJ(B,C,fill,BC);

 34:   MatProductSetAlgorithm(D,"sorted"); /* set alg for D = A*BC */
 35:   MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,BC,fill,D);
 36:   D->product->alg = alg; /* resume original algorithm for D */

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

 41:   matmatmatmult->BC      = BC;
 42:   matmatmatmult->destroy = D->ops->destroy;
 43:   d                      = (Mat_SeqAIJ*)D->data;
 44:   d->matmatmatmult       = matmatmatmult;

 46:   D->ops->matmatmultnumeric = MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ;
 47:   D->ops->destroy           = MatDestroy_SeqAIJ_MatMatMatMult;
 48:   return(0);
 49: }

 51: PetscErrorCode MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,Mat D)
 52: {
 53:   PetscErrorCode    ierr;
 54:   Mat_SeqAIJ        *d            =(Mat_SeqAIJ*)D->data;
 55:   Mat_MatMatMatMult *matmatmatmult=d->matmatmatmult;
 56:   Mat               BC            = matmatmatmult->BC;

 59:   (BC->ops->matmultnumeric)(B,C,BC);
 60:   (D->ops->matmultnumeric)(A,BC,D);
 61:   return(0);
 62: }