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;
20: PetscInitialize(&argc, &argv, NULL, help);
21: comm = PETSC_COMM_WORLD;
22: MPI_Comm_rank(comm, &rank);
23: MPI_Comm_size(comm, &size);
24: M = 10;
25: N = 10;
26: Z = 10;
27: dof = 10;
29: PetscOptionsGetBool(NULL, NULL, "-test_3D", &Test_3D, NULL);
30: PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL);
31: PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL);
32: PetscOptionsGetInt(NULL, NULL, "-Z", &Z, NULL);
33: /* Set up distributed array for fine grid */
34: if (!Test_3D) {
35: 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: 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: DMSetFromOptions(coarsedm);
40: DMSetUp(coarsedm);
42: /* This makes sure the coarse DMDA has the same partition as the fine DMDA */
43: DMRefine(coarsedm, PetscObjectComm((PetscObject)coarsedm), &finedm);
45: /*------------------------------------------------------------*/
46: DMSetMatType(finedm, MATAIJ);
47: DMCreateMatrix(finedm, &A);
49: /* set val=one to A */
50: if (size == 1) {
51: MatGetRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg);
52: if (flg) {
53: MatSeqAIJGetArray(A, &array);
54: for (i = 0; i < ia[nrows]; i++) array[i] = one;
55: MatSeqAIJRestoreArray(A, &array);
56: }
57: MatRestoreRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg);
58: } else {
59: Mat AA, AB;
60: MatMPIAIJGetSeqAIJ(A, &AA, &AB, NULL);
61: MatGetRowIJ(AA, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg);
62: if (flg) {
63: MatSeqAIJGetArray(AA, &array);
64: for (i = 0; i < ia[nrows]; i++) array[i] = one;
65: MatSeqAIJRestoreArray(AA, &array);
66: }
67: MatRestoreRowIJ(AA, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg);
68: MatGetRowIJ(AB, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg);
69: if (flg) {
70: MatSeqAIJGetArray(AB, &array);
71: for (i = 0; i < ia[nrows]; i++) array[i] = one;
72: MatSeqAIJRestoreArray(AB, &array);
73: }
74: MatRestoreRowIJ(AB, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg);
75: }
76: /* Create interpolation between the fine and coarse grids */
77: DMCreateInterpolation(coarsedm, finedm, &P, NULL);
79: /* Test P^T * A * P - MatPtAP() */
80: /*------------------------------*/
81: /* (1) Developer API */
82: MatProductCreate(A, P, NULL, &C);
83: MatProductSetType(C, MATPRODUCT_PtAP);
84: MatProductSetAlgorithm(C, "allatonce");
85: MatProductSetFill(C, PETSC_DEFAULT);
86: MatProductSetFromOptions(C);
87: MatProductSymbolic(C);
88: MatProductNumeric(C);
89: MatProductNumeric(C); /* Test reuse of symbolic C */
91: { /* Test MatProductView() */
92: PetscViewer viewer;
93: PetscViewerASCIIOpen(comm, NULL, &viewer);
94: PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO);
95: MatProductView(C, viewer);
96: PetscViewerPopFormat(viewer);
97: PetscViewerDestroy(&viewer);
98: }
100: MatPtAPMultEqual(A, P, C, 10, &flg);
102: MatDestroy(&C);
104: /* (2) User API */
105: 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: MatScale(A, alpha);
111: MatPtAP(A, P, MAT_REUSE_MATRIX, fill, &C);
112: }
114: /* Free intermediate data structures created for reuse of C=Pt*A*P */
115: MatProductClear(C);
117: MatPtAPMultEqual(A, P, C, 10, &flg);
120: MatDestroy(&C);
121: MatDestroy(&A);
122: MatDestroy(&P);
123: DMDestroy(&finedm);
124: DMDestroy(&coarsedm);
125: 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*/