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;
     10:   PetscErrorCode    ierr;

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

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

     26:   MatCheckProduct(D,5);
     27:   if (D->product->data) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Product data not empty");
     28:   MatCreate(PETSC_COMM_SELF,&BC);
     29:   MatMatMultSymbolic_SeqAIJ_SeqAIJ(B,C,fill,BC);

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

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

     44:   D->ops->matmatmultnumeric = MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ;
     45:   return(0);
     46: }

     48: PetscErrorCode MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,Mat D)
     49: {
     50:   PetscErrorCode    ierr;
     51:   Mat_MatMatMatMult *matmatmatmult;
     52:   Mat               BC;

     55:   MatCheckProduct(D,4);
     56:   if (!D->product->data) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Product data empty");
     57:   matmatmatmult = (Mat_MatMatMatMult*)D->product->data;
     58:   BC = matmatmatmult->BC;
     59:   if (!BC) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Missing BC mat");
     60:   if (!BC->ops->matmultnumeric) SETERRQ(PetscObjectComm((PetscObject)BC),PETSC_ERR_PLIB,"Missing numeric operation");
     61:   (*BC->ops->matmultnumeric)(B,C,BC);
     62:   if (!D->ops->matmultnumeric) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Missing numeric operation");
     63:   (*D->ops->matmultnumeric)(A,BC,D);
     64:   return(0);
     65: }