1 static char help[] = "Tests MatPtAP() for MPIMAIJ and MPIAIJ \n ";
2
3 #include <petscdmda.h>
4
main(int argc,char ** argv)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;
18
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;
28
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));
41
42 /* This makes sure the coarse DMDA has the same partition as the fine DMDA */
43 PetscCall(DMRefine(coarsedm, PetscObjectComm((PetscObject)coarsedm), &finedm));
44
45 /*------------------------------------------------------------*/
46 PetscCall(DMSetMatType(finedm, MATAIJ));
47 PetscCall(DMCreateMatrix(finedm, &A));
48
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));
78
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 */
90
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 }
99
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));
103
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 }
113
114 /* Free intermediate data structures created for reuse of C=Pt*A*P */
115 PetscCall(MatProductClear(C));
116
117 PetscCall(MatPtAPMultEqual(A, P, C, 10, &flg));
118 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in MatPtAP");
119
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 }
128
129 /*TEST
130
131 test:
132 args: -M 10 -N 10 -Z 10
133 output_file: output/empty.out
134
135 test:
136 suffix: allatonce
137 nsize: 4
138 args: -M 10 -N 10 -Z 10
139 output_file: output/ex89_2.out
140
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
146
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
152
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
158
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
164
165 TEST*/
166