Actual source code: matmatmatmult.c

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

 11:   MatDestroy(&matmatmatmult->BC);
 12:   PetscFree(matmatmatmult);
 13:   return 0;
 14: }

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

 22:   MatCheckProduct(D, 5);
 24:   MatCreate(PETSC_COMM_SELF, &BC);
 25:   MatMatMultSymbolic_SeqAIJ_SeqAIJ(B, C, fill, BC);

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

 33:   /* create struct Mat_MatMatMatMult and attached it to D */
 35:   PetscNew(&matmatmatmult);
 36:   matmatmatmult->BC   = BC;
 37:   D->product->data    = matmatmatmult;
 38:   D->product->destroy = MatDestroy_SeqAIJ_MatMatMatMult;

 40:   D->ops->matmatmultnumeric = MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ;
 41:   return 0;
 42: }

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

 49:   MatCheckProduct(D, 4);
 51:   matmatmatmult = (Mat_MatMatMatMult *)D->product->data;
 52:   BC            = matmatmatmult->BC;
 54:   (*BC->ops->matmultnumeric)(B, C, BC);
 55:   (*D->ops->matmultnumeric)(A, BC, D);
 56:   return 0;
 57: }