Actual source code: ex89.c
1: static char help[] = "Tests MatPtAP() for MPIMAIJ and MPIAIJ \n ";
3: #include <petscdmda.h>
5: int main(int argc, char **argv)
6: {
7: DM coarsedm, finedm;
8: PetscMPIInt size, rank;
9: PetscInt M, N, Z, i, nrows;
10: PetscScalar one = 1.0;
11: PetscReal fill = 2.0;
12: Mat A, P, C;
13: PetscScalar *array, alpha;
14: PetscBool Test_3D = PETSC_FALSE, flg;
15: const PetscInt *ia, *ja;
16: PetscInt dof;
17: MPI_Comm comm;
19: PetscFunctionBeginUser;
20: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
21: comm = PETSC_COMM_WORLD;
22: PetscCallMPI(MPI_Comm_rank(comm, &rank));
23: PetscCallMPI(MPI_Comm_size(comm, &size));
24: M = 10;
25: N = 10;
26: Z = 10;
27: dof = 10;
29: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_3D", &Test_3D, NULL));
30: PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
31: PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL));
32: PetscCall(PetscOptionsGetInt(NULL, NULL, "-Z", &Z, NULL));
33: /* Set up distributed array for fine grid */
34: if (!Test_3D) {
35: PetscCall(DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, M, N, PETSC_DECIDE, PETSC_DECIDE, dof, 1, NULL, NULL, &coarsedm));
36: } else {
37: PetscCall(DMDACreate3d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, M, N, Z, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, dof, 1, NULL, NULL, NULL, &coarsedm));
38: }
39: PetscCall(DMSetFromOptions(coarsedm));
40: PetscCall(DMSetUp(coarsedm));
42: /* This makes sure the coarse DMDA has the same partition as the fine DMDA */
43: PetscCall(DMRefine(coarsedm, PetscObjectComm((PetscObject)coarsedm), &finedm));
45: /*------------------------------------------------------------*/
46: PetscCall(DMSetMatType(finedm, MATAIJ));
47: PetscCall(DMCreateMatrix(finedm, &A));
49: /* set val=one to A */
50: if (size == 1) {
51: PetscCall(MatGetRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg));
52: if (flg) {
53: PetscCall(MatSeqAIJGetArray(A, &array));
54: for (i = 0; i < ia[nrows]; i++) array[i] = one;
55: PetscCall(MatSeqAIJRestoreArray(A, &array));
56: }
57: PetscCall(MatRestoreRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg));
58: } else {
59: Mat AA, AB;
60: PetscCall(MatMPIAIJGetSeqAIJ(A, &AA, &AB, NULL));
61: PetscCall(MatGetRowIJ(AA, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg));
62: if (flg) {
63: PetscCall(MatSeqAIJGetArray(AA, &array));
64: for (i = 0; i < ia[nrows]; i++) array[i] = one;
65: PetscCall(MatSeqAIJRestoreArray(AA, &array));
66: }
67: PetscCall(MatRestoreRowIJ(AA, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg));
68: PetscCall(MatGetRowIJ(AB, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg));
69: if (flg) {
70: PetscCall(MatSeqAIJGetArray(AB, &array));
71: for (i = 0; i < ia[nrows]; i++) array[i] = one;
72: PetscCall(MatSeqAIJRestoreArray(AB, &array));
73: }
74: PetscCall(MatRestoreRowIJ(AB, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg));
75: }
76: /* Create interpolation between the fine and coarse grids */
77: PetscCall(DMCreateInterpolation(coarsedm, finedm, &P, NULL));
79: /* Test P^T * A * P - MatPtAP() */
80: /*------------------------------*/
81: /* (1) Developer API */
82: PetscCall(MatProductCreate(A, P, NULL, &C));
83: PetscCall(MatProductSetType(C, MATPRODUCT_PtAP));
84: PetscCall(MatProductSetAlgorithm(C, "allatonce"));
85: PetscCall(MatProductSetFill(C, PETSC_DEFAULT));
86: PetscCall(MatProductSetFromOptions(C));
87: PetscCall(MatProductSymbolic(C));
88: PetscCall(MatProductNumeric(C));
89: PetscCall(MatProductNumeric(C)); /* Test reuse of symbolic C */
91: { /* Test MatProductView() */
92: PetscViewer viewer;
93: PetscCall(PetscViewerASCIIOpen(comm, NULL, &viewer));
94: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO));
95: PetscCall(MatProductView(C, viewer));
96: PetscCall(PetscViewerPopFormat(viewer));
97: PetscCall(PetscViewerDestroy(&viewer));
98: }
100: PetscCall(MatPtAPMultEqual(A, P, C, 10, &flg));
101: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatProduct_PtAP");
102: PetscCall(MatDestroy(&C));
104: /* (2) User API */
105: PetscCall(MatPtAP(A, P, MAT_INITIAL_MATRIX, fill, &C));
106: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
107: alpha = 1.0;
108: for (i = 0; i < 1; i++) {
109: alpha -= 0.1;
110: PetscCall(MatScale(A, alpha));
111: PetscCall(MatPtAP(A, P, MAT_REUSE_MATRIX, fill, &C));
112: }
114: /* Free intermediate data structures created for reuse of C=Pt*A*P */
115: PetscCall(MatProductClear(C));
117: PetscCall(MatPtAPMultEqual(A, P, C, 10, &flg));
118: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatPtAP");
120: PetscCall(MatDestroy(&C));
121: PetscCall(MatDestroy(&A));
122: PetscCall(MatDestroy(&P));
123: PetscCall(DMDestroy(&finedm));
124: PetscCall(DMDestroy(&coarsedm));
125: PetscCall(PetscFinalize());
126: return 0;
127: }
129: /*TEST
131: test:
132: args: -M 10 -N 10 -Z 10
133: output_file: output/ex89_1.out
135: test:
136: suffix: allatonce
137: nsize: 4
138: args: -M 10 -N 10 -Z 10
139: output_file: output/ex89_2.out
141: test:
142: suffix: allatonce_merged
143: nsize: 4
144: args: -M 10 -M 5 -M 10 -mat_product_algorithm allatonce_merged
145: output_file: output/ex89_3.out
147: test:
148: suffix: nonscalable_3D
149: nsize: 4
150: args: -M 10 -M 5 -M 10 -test_3D 1 -mat_product_algorithm nonscalable
151: output_file: output/ex89_4.out
153: test:
154: suffix: allatonce_merged_3D
155: nsize: 4
156: args: -M 10 -M 5 -M 10 -test_3D 1 -mat_product_algorithm allatonce_merged
157: output_file: output/ex89_3.out
159: test:
160: suffix: nonscalable
161: nsize: 4
162: args: -M 10 -N 10 -Z 10 -mat_product_algorithm nonscalable
163: output_file: output/ex89_5.out
165: TEST*/