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: 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: }