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