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