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> /*I "petscmat.h" I*/
6
MatProductCtxDestroy_SeqAIJ_MatMatMatMult(PetscCtxRt data)7 static PetscErrorCode MatProductCtxDestroy_SeqAIJ_MatMatMatMult(PetscCtxRt data)
8 {
9 MatProductCtx_MatMatMatMult *matmatmatmult = *(MatProductCtx_MatMatMatMult **)data;
10
11 PetscFunctionBegin;
12 PetscCall(MatDestroy(&matmatmatmult->BC));
13 PetscCall(PetscFree(matmatmatmult));
14 PetscFunctionReturn(PETSC_SUCCESS);
15 }
16
MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,PetscReal fill,Mat D)17 PetscErrorCode MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(Mat A, Mat B, Mat C, PetscReal fill, Mat D)
18 {
19 Mat BC;
20 MatProductCtx_MatMatMatMult *matmatmatmult;
21 char *alg;
22
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));
28
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));
34
35 /* create struct MatProductCtx_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 = MatProductCtxDestroy_SeqAIJ_MatMatMatMult;
41
42 D->ops->matmatmultnumeric = MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ;
43 PetscFunctionReturn(PETSC_SUCCESS);
44 }
45
MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C,Mat D)46 PetscErrorCode MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(Mat A, Mat B, Mat C, Mat D)
47 {
48 MatProductCtx_MatMatMatMult *matmatmatmult;
49 Mat BC;
50
51 PetscFunctionBegin;
52 MatCheckProduct(D, 4);
53 PetscCheck(D->product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data empty");
54 matmatmatmult = (MatProductCtx_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 }
61