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: static PetscErrorCode MatDestroy_SeqAIJ_MatMatMatMult(void *data)
  8: {
  9:   Mat_MatMatMatMult *matmatmatmult = (Mat_MatMatMatMult *)data;

 11:   PetscFunctionBegin;
 12:   PetscCall(MatDestroy(&matmatmatmult->BC));
 13:   PetscCall(PetscFree(matmatmatmult));
 14:   PetscFunctionReturn(PETSC_SUCCESS);
 15: }

 17: PetscErrorCode MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(Mat A, Mat B, Mat C, PetscReal fill, Mat D)
 18: {
 19:   Mat                BC;
 20:   Mat_MatMatMatMult *matmatmatmult;
 21:   char              *alg;

 23:   PetscFunctionBegin;
 24:   MatCheckProduct(D, 5);
 25:   PetscCheck(!D->product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data not empty");
 26:   PetscCall(MatCreate(PETSC_COMM_SELF, &BC));
 27:   PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ(B, C, fill, BC));

 29:   PetscCall(PetscStrallocpy(D->product->alg, &alg));
 30:   PetscCall(MatProductSetAlgorithm(D, "sorted")); /* set alg for D = A*BC */
 31:   PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ(A, BC, fill, D));
 32:   PetscCall(MatProductSetAlgorithm(D, alg)); /* resume original algorithm */
 33:   PetscCall(PetscFree(alg));

 35:   /* create struct Mat_MatMatMatMult and attached it to D */
 36:   PetscCheck(!D->product->data, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Not yet coded");
 37:   PetscCall(PetscNew(&matmatmatmult));
 38:   matmatmatmult->BC   = BC;
 39:   D->product->data    = matmatmatmult;
 40:   D->product->destroy = MatDestroy_SeqAIJ_MatMatMatMult;

 42:   D->ops->matmatmultnumeric = MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ;
 43:   PetscFunctionReturn(PETSC_SUCCESS);
 44: }

 46: PetscErrorCode MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(Mat A, Mat B, Mat C, Mat D)
 47: {
 48:   Mat_MatMatMatMult *matmatmatmult;
 49:   Mat                BC;

 51:   PetscFunctionBegin;
 52:   MatCheckProduct(D, 4);
 53:   PetscCheck(D->product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data empty");
 54:   matmatmatmult = (Mat_MatMatMatMult *)D->product->data;
 55:   BC            = matmatmatmult->BC;
 56:   PetscCheck(BC, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Missing BC mat");
 57:   PetscCall((*BC->ops->matmultnumeric)(B, C, BC));
 58:   PetscCall((*D->ops->matmultnumeric)(A, BC, D));
 59:   PetscFunctionReturn(PETSC_SUCCESS);
 60: }