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*/