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;
55: (*BC->ops->matmultnumeric)(B,C,BC);
57: (*D->ops->matmultnumeric)(A,BC,D);
58: return 0;
59: }