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 %" PetscInt_FMT "\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*/