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;
34: 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: 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: DMSetFromOptions(da);
42: 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: 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: PetscLogStageRegister("Assembly", &stage);
63: PetscLogStagePush(stage);
64: 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: 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: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
108: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
109: PetscLogStagePop();
111: /*
112: Create parallel vectors compatible with the DMDA.
113: */
114: DMCreateGlobalVector(da, &u);
115: VecDuplicate(u, &b);
116: 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: PetscOptionsGetBool(NULL, NULL, "-random_exact_sol", &flg, NULL);
125: if (flg) {
126: PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
127: PetscRandomSetFromOptions(rctx);
128: VecSetRandom(u, rctx);
129: PetscRandomDestroy(&rctx);
130: } else {
131: VecSet(u, 1.);
132: }
133: MatMult(A, u, b);
135: /*
136: View the exact solution vector if desired
137: */
138: flg = PETSC_FALSE;
139: PetscOptionsGetBool(NULL, NULL, "-view_exact_sol", &flg, NULL);
140: if (flg) 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: 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: 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: KSPSetFromOptions(ksp);
166: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167: Solve the linear system
168: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
170: KSPSolve(ksp, b, x);
172: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
173: Check solution and clean up
174: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
176: /*
177: Check the error
178: */
179: VecAXPY(x, -1., u);
180: VecNorm(x, NORM_2, &norm);
181: 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: 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: KSPDestroy(&ksp);
195: VecDestroy(&u);
196: VecDestroy(&x);
197: VecDestroy(&b);
198: MatDestroy(&A);
199: 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: 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*/