Actual source code: ex46.c

      1: static char help[] = "Solves a linear system in parallel with KSP and DM.\n\
      2: Compare this to ex2 which solves the same problem without a DM.\n\n";

      4: /*
      5:   Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
      6:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
      7:   automatically includes:
      8:      petscsys.h       - base PETSc routines   petscvec.h - vectors
      9:      petscmat.h - matrices
     10:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
     11:      petscviewer.h - viewers               petscpc.h  - preconditioners
     12: */
     13: #include <petscdm.h>
     14: #include <petscdmda.h>
     15: #include <petscksp.h>

     17: int main(int argc, char **argv)
     18: {
     19:   DM            da;      /* distributed array */
     20:   Vec           x, b, u; /* approx solution, RHS, exact solution */
     21:   Mat           A;       /* linear system matrix */
     22:   KSP           ksp;     /* linear solver context */
     23:   PetscRandom   rctx;    /* random number generator context */
     24:   PetscReal     norm;    /* norm of solution error */
     25:   PetscInt      i, j, its;
     26:   PetscBool     flg = PETSC_FALSE;
     27:   PetscLogStage stage;
     28:   DMDALocalInfo info;

     30:   PetscFunctionBeginUser;
     31:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
     32:   /*
     33:      Create distributed array to handle parallel distribution.
     34:      The problem size will default to 8 by 7, but this can be
     35:      changed using -da_grid_x M -da_grid_y N
     36:   */
     37:   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));
     38:   PetscCall(DMSetFromOptions(da));
     39:   PetscCall(DMSetUp(da));

     41:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     42:          Compute the matrix and right-hand-side vector that define
     43:          the linear system, Ax = b.
     44:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     45:   /*
     46:      Create parallel matrix preallocated according to the DMDA, format AIJ by default.
     47:      To use symmetric storage, run with -dm_mat_type sbaij -mat_ignore_lower_triangular
     48:   */
     49:   PetscCall(DMCreateMatrix(da, &A));

     51:   /*
     52:      Set matrix elements for the 2-D, five-point stencil in parallel.
     53:       - Each processor needs to insert only elements that it owns
     54:         locally (but any non-local elements will be sent to the
     55:         appropriate processor during matrix assembly).
     56:       - Rows and columns are specified by the stencil
     57:       - Entries are normalized for a domain [0,1]x[0,1]
     58:    */
     59:   PetscCall(PetscLogStageRegister("Assembly", &stage));
     60:   PetscCall(PetscLogStagePush(stage));
     61:   PetscCall(DMDAGetLocalInfo(da, &info));
     62:   for (j = info.ys; j < info.ys + info.ym; j++) {
     63:     for (i = info.xs; i < info.xs + info.xm; i++) {
     64:       PetscReal   hx = 1. / info.mx, hy = 1. / info.my;
     65:       MatStencil  row = {0}, col[5] = {{0}};
     66:       PetscScalar v[5];
     67:       PetscInt    ncols = 0;
     68:       row.j             = j;
     69:       row.i             = i;
     70:       col[ncols].j      = j;
     71:       col[ncols].i      = i;
     72:       v[ncols++]        = 2 * (hx / hy + hy / hx);
     73:       /* boundaries */
     74:       if (i > 0) {
     75:         col[ncols].j = j;
     76:         col[ncols].i = i - 1;
     77:         v[ncols++]   = -hy / hx;
     78:       }
     79:       if (i < info.mx - 1) {
     80:         col[ncols].j = j;
     81:         col[ncols].i = i + 1;
     82:         v[ncols++]   = -hy / hx;
     83:       }
     84:       if (j > 0) {
     85:         col[ncols].j = j - 1;
     86:         col[ncols].i = i;
     87:         v[ncols++]   = -hx / hy;
     88:       }
     89:       if (j < info.my - 1) {
     90:         col[ncols].j = j + 1;
     91:         col[ncols].i = i;
     92:         v[ncols++]   = -hx / hy;
     93:       }
     94:       PetscCall(MatSetValuesStencil(A, 1, &row, ncols, col, v, INSERT_VALUES));
     95:     }
     96:   }

     98:   /*
     99:      Assemble matrix, using the 2-step process:
    100:        MatAssemblyBegin(), MatAssemblyEnd()
    101:      Computations can be done while messages are in transition
    102:      by placing code between these two statements.
    103:   */
    104:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    105:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    106:   PetscCall(PetscLogStagePop());

    108:   /*
    109:      Create parallel vectors compatible with the DMDA.
    110:   */
    111:   PetscCall(DMCreateGlobalVector(da, &u));
    112:   PetscCall(VecDuplicate(u, &b));
    113:   PetscCall(VecDuplicate(u, &x));

    115:   /*
    116:      Set exact solution; then compute right-hand-side vector.
    117:      By default we use an exact solution of a vector with all
    118:      elements of 1.0;  Alternatively, using the runtime option
    119:      -random_sol forms a solution vector with random components.
    120:   */
    121:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-random_exact_sol", &flg, NULL));
    122:   if (flg) {
    123:     PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
    124:     PetscCall(PetscRandomSetFromOptions(rctx));
    125:     PetscCall(VecSetRandom(u, rctx));
    126:     PetscCall(PetscRandomDestroy(&rctx));
    127:   } else {
    128:     PetscCall(VecSet(u, 1.));
    129:   }
    130:   PetscCall(MatMult(A, u, b));

    132:   /*
    133:      View the exact solution vector if desired
    134:   */
    135:   flg = PETSC_FALSE;
    136:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_exact_sol", &flg, NULL));
    137:   if (flg) PetscCall(VecView(u, PETSC_VIEWER_STDOUT_WORLD));

    139:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    140:                 Create the linear solver and set various options
    141:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    143:   /*
    144:      Create linear solver context
    145:   */
    146:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));

    148:   /*
    149:      Set operators. Here the matrix that defines the linear system
    150:      also serves as the matrix from which the preconditioner is constructed.
    151:   */
    152:   PetscCall(KSPSetOperators(ksp, A, A));

    154:   /*
    155:     Set runtime options, e.g.,
    156:         -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
    157:     These options will override those specified above as long as
    158:     KSPSetFromOptions() is called _after_ any other customization
    159:     routines.
    160:   */
    161:   PetscCall(KSPSetFromOptions(ksp));

    163:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    164:                       Solve the linear system
    165:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    167:   PetscCall(KSPSolve(ksp, b, x));

    169:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    170:                       Check solution and clean up
    171:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    173:   /*
    174:      Check the error
    175:   */
    176:   PetscCall(VecAXPY(x, -1., u));
    177:   PetscCall(VecNorm(x, NORM_2, &norm));
    178:   PetscCall(KSPGetIterationNumber(ksp, &its));

    180:   /*
    181:      Print convergence information.  PetscPrintf() produces a single
    182:      print statement from all processes that share a communicator.
    183:      An alternative is PetscFPrintf(), which prints to a file.
    184:   */
    185:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g iterations %" PetscInt_FMT "\n", (double)norm, its));

    187:   /*
    188:      Free work space.  All PETSc objects should be destroyed when they
    189:      are no longer needed.
    190:   */
    191:   PetscCall(KSPDestroy(&ksp));
    192:   PetscCall(VecDestroy(&u));
    193:   PetscCall(VecDestroy(&x));
    194:   PetscCall(VecDestroy(&b));
    195:   PetscCall(MatDestroy(&A));
    196:   PetscCall(DMDestroy(&da));

    198:   /*
    199:      Always call PetscFinalize() before exiting a program.  This routine
    200:        - finalizes the PETSc libraries as well as MPI
    201:        - provides summary and diagnostic information if certain runtime
    202:          options are chosen (e.g., -log_view).
    203:   */
    204:   PetscCall(PetscFinalize());
    205:   return 0;
    206: }

    208: /*TEST

    210:    testset:
    211:       requires: cuda
    212:       args: -dm_mat_type aijcusparse -dm_vec_type cuda -random_exact_sol
    213:       output_file: output/ex46_aijcusparse.out

    215:       test:
    216:         suffix: aijcusparse
    217:         args: -use_gpu_aware_mpi 0
    218:       test:
    219:         suffix: aijcusparse_mpi_gpu_aware

    221:    testset:
    222:       requires: cuda
    223:       args: -dm_mat_type aijcusparse -dm_vec_type cuda -random_exact_sol -pc_type ilu -pc_factor_mat_solver_type cusparse
    224:       output_file: output/ex46_aijcusparse_2.out
    225:       test:
    226:         suffix: aijcusparse_2
    227:         args: -use_gpu_aware_mpi 0
    228:       test:
    229:         suffix: aijcusparse_2_mpi_gpu_aware

    231: TEST*/