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