Actual source code: ex46.c
2: static char help[] = "Solves a linear system in parallel with KSP and DM.\n\
3: Compare this to ex2 which solves the same problem without a DM.\n\n";
5: /*
6: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
7: Include "petscksp.h" so that we can use KSP solvers. Note that this file
8: automatically includes:
9: petscsys.h - base PETSc routines petscvec.h - vectors
10: petscmat.h - matrices
11: petscis.h - index sets petscksp.h - Krylov subspace methods
12: petscviewer.h - viewers petscpc.h - preconditioners
13: */
14: #include <petscdm.h>
15: #include <petscdmda.h>
16: #include <petscksp.h>
18: int main(int argc, char **argv)
19: {
20: DM da; /* distributed array */
21: Vec x, b, u; /* approx solution, RHS, exact solution */
22: Mat A; /* linear system matrix */
23: KSP ksp; /* linear solver context */
24: PetscRandom rctx; /* random number generator context */
25: PetscReal norm; /* norm of solution error */
26: PetscInt i, j, its;
27: PetscBool flg = PETSC_FALSE;
28: #if defined(PETSC_USE_LOG)
29: PetscLogStage stage;
30: #endif
31: DMDALocalInfo info;
33: PetscFunctionBeginUser;
34: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
35: /*
36: Create distributed array to handle parallel distribution.
37: The problem size will default to 8 by 7, but this can be
38: changed using -da_grid_x M -da_grid_y N
39: */
40: PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 8, 7, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da));
41: PetscCall(DMSetFromOptions(da));
42: PetscCall(DMSetUp(da));
44: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
45: Compute the matrix and right-hand-side vector that define
46: the linear system, Ax = b.
47: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
48: /*
49: Create parallel matrix preallocated according to the DMDA, format AIJ by default.
50: To use symmetric storage, run with -dm_mat_type sbaij -mat_ignore_lower_triangular
51: */
52: PetscCall(DMCreateMatrix(da, &A));
54: /*
55: Set matrix elements for the 2-D, five-point stencil in parallel.
56: - Each processor needs to insert only elements that it owns
57: locally (but any non-local elements will be sent to the
58: appropriate processor during matrix assembly).
59: - Rows and columns are specified by the stencil
60: - Entries are normalized for a domain [0,1]x[0,1]
61: */
62: PetscCall(PetscLogStageRegister("Assembly", &stage));
63: PetscCall(PetscLogStagePush(stage));
64: PetscCall(DMDAGetLocalInfo(da, &info));
65: for (j = info.ys; j < info.ys + info.ym; j++) {
66: for (i = info.xs; i < info.xs + info.xm; i++) {
67: PetscReal hx = 1. / info.mx, hy = 1. / info.my;
68: MatStencil row = {0}, col[5] = {{0}};
69: PetscScalar v[5];
70: PetscInt ncols = 0;
71: row.j = j;
72: row.i = i;
73: col[ncols].j = j;
74: col[ncols].i = i;
75: v[ncols++] = 2 * (hx / hy + hy / hx);
76: /* boundaries */
77: if (i > 0) {
78: col[ncols].j = j;
79: col[ncols].i = i - 1;
80: v[ncols++] = -hy / hx;
81: }
82: if (i < info.mx - 1) {
83: col[ncols].j = j;
84: col[ncols].i = i + 1;
85: v[ncols++] = -hy / hx;
86: }
87: if (j > 0) {
88: col[ncols].j = j - 1;
89: col[ncols].i = i;
90: v[ncols++] = -hx / hy;
91: }
92: if (j < info.my - 1) {
93: col[ncols].j = j + 1;
94: col[ncols].i = i;
95: v[ncols++] = -hx / hy;
96: }
97: PetscCall(MatSetValuesStencil(A, 1, &row, ncols, col, v, INSERT_VALUES));
98: }
99: }
101: /*
102: Assemble matrix, using the 2-step process:
103: MatAssemblyBegin(), MatAssemblyEnd()
104: Computations can be done while messages are in transition
105: by placing code between these two statements.
106: */
107: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
108: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
109: PetscCall(PetscLogStagePop());
111: /*
112: Create parallel vectors compatible with the DMDA.
113: */
114: PetscCall(DMCreateGlobalVector(da, &u));
115: PetscCall(VecDuplicate(u, &b));
116: PetscCall(VecDuplicate(u, &x));
118: /*
119: Set exact solution; then compute right-hand-side vector.
120: By default we use an exact solution of a vector with all
121: elements of 1.0; Alternatively, using the runtime option
122: -random_sol forms a solution vector with random components.
123: */
124: PetscCall(PetscOptionsGetBool(NULL, NULL, "-random_exact_sol", &flg, NULL));
125: if (flg) {
126: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
127: PetscCall(PetscRandomSetFromOptions(rctx));
128: PetscCall(VecSetRandom(u, rctx));
129: PetscCall(PetscRandomDestroy(&rctx));
130: } else {
131: PetscCall(VecSet(u, 1.));
132: }
133: PetscCall(MatMult(A, u, b));
135: /*
136: View the exact solution vector if desired
137: */
138: flg = PETSC_FALSE;
139: PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_exact_sol", &flg, NULL));
140: if (flg) PetscCall(VecView(u, PETSC_VIEWER_STDOUT_WORLD));
142: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143: Create the linear solver and set various options
144: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
146: /*
147: Create linear solver context
148: */
149: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
151: /*
152: Set operators. Here the matrix that defines the linear system
153: also serves as the preconditioning matrix.
154: */
155: PetscCall(KSPSetOperators(ksp, A, A));
157: /*
158: Set runtime options, e.g.,
159: -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
160: These options will override those specified above as long as
161: KSPSetFromOptions() is called _after_ any other customization
162: routines.
163: */
164: PetscCall(KSPSetFromOptions(ksp));
166: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167: Solve the linear system
168: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
170: PetscCall(KSPSolve(ksp, b, x));
172: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
173: Check solution and clean up
174: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
176: /*
177: Check the error
178: */
179: PetscCall(VecAXPY(x, -1., u));
180: PetscCall(VecNorm(x, NORM_2, &norm));
181: PetscCall(KSPGetIterationNumber(ksp, &its));
183: /*
184: Print convergence information. PetscPrintf() produces a single
185: print statement from all processes that share a communicator.
186: An alternative is PetscFPrintf(), which prints to a file.
187: */
188: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g iterations %" PetscInt_FMT "\n", (double)norm, its));
190: /*
191: Free work space. All PETSc objects should be destroyed when they
192: are no longer needed.
193: */
194: PetscCall(KSPDestroy(&ksp));
195: PetscCall(VecDestroy(&u));
196: PetscCall(VecDestroy(&x));
197: PetscCall(VecDestroy(&b));
198: PetscCall(MatDestroy(&A));
199: PetscCall(DMDestroy(&da));
201: /*
202: Always call PetscFinalize() before exiting a program. This routine
203: - finalizes the PETSc libraries as well as MPI
204: - provides summary and diagnostic information if certain runtime
205: options are chosen (e.g., -log_view).
206: */
207: PetscCall(PetscFinalize());
208: return 0;
209: }
211: /*TEST
213: testset:
214: requires: cuda
215: args: -dm_mat_type aijcusparse -dm_vec_type cuda -random_exact_sol
216: output_file: output/ex46_aijcusparse.out
218: test:
219: suffix: aijcusparse
220: args: -use_gpu_aware_mpi 0
221: test:
222: suffix: aijcusparse_mpi_gpu_aware
224: testset:
225: requires: cuda
226: args: -dm_mat_type aijcusparse -dm_vec_type cuda -random_exact_sol -pc_type ilu -pc_factor_mat_solver_type cusparse
227: output_file: output/ex46_aijcusparse_2.out
228: test:
229: suffix: aijcusparse_2
230: args: -use_gpu_aware_mpi 0
231: test:
232: suffix: aijcusparse_2_mpi_gpu_aware
234: TEST*/