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