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