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