Actual source code: ex2.c


  2: static char help[] = "Solves a linear system in parallel with KSP.\n\
  3: Input parameters include:\n\
  4:   -view_exact_sol   : write exact solution vector to stdout\n\
  5:   -m <mesh_x>       : number of mesh points in x-direction\n\
  6:   -n <mesh_y>       : number of mesh points in y-direction\n\n";

  8: /*
  9:   Include "petscksp.h" so that we can use KSP solvers.
 10: */
 11: #include <petscksp.h>

 13: int main(int argc,char **args)
 14: {
 15:   Vec            x,b,u;    /* approx solution, RHS, exact solution */
 16:   Mat            A;        /* linear system matrix */
 17:   KSP            ksp;      /* linear solver context */
 18:   PetscReal      norm;     /* norm of solution error */
 19:   PetscInt       i,j,Ii,J,Istart,Iend,m = 8,n = 7,its;
 20:   PetscBool      flg;
 21:   PetscScalar    v;

 23:   PetscInitialize(&argc,&args,(char*)0,help);
 24:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
 25:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 26:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 27:          Compute the matrix and right-hand-side vector that define
 28:          the linear system, Ax = b.
 29:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 30:   /*
 31:      Create parallel matrix, specifying only its global dimensions.
 32:      When using MatCreate(), the matrix format can be specified at
 33:      runtime. Also, the parallel partitioning of the matrix is
 34:      determined by PETSc at runtime.

 36:      Performance tuning note:  For problems of substantial size,
 37:      preallocation of matrix memory is crucial for attaining good
 38:      performance. See the matrix chapter of the users manual for details.
 39:   */
 40:   MatCreate(PETSC_COMM_WORLD,&A);
 41:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
 42:   MatSetFromOptions(A);
 43:   MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);
 44:   MatSeqAIJSetPreallocation(A,5,NULL);
 45:   MatSeqSBAIJSetPreallocation(A,1,5,NULL);
 46:   MatMPISBAIJSetPreallocation(A,1,5,NULL,5,NULL);
 47:   MatMPISELLSetPreallocation(A,5,NULL,5,NULL);
 48:   MatSeqSELLSetPreallocation(A,5,NULL);

 50:   /*
 51:      Currently, all PETSc parallel matrix formats are partitioned by
 52:      contiguous chunks of rows across the processors.  Determine which
 53:      rows of the matrix are locally owned.
 54:   */
 55:   MatGetOwnershipRange(A,&Istart,&Iend);

 57:   /*
 58:      Set matrix elements for the 2-D, five-point stencil in parallel.
 59:       - Each processor needs to insert only elements that it owns
 60:         locally (but any non-local elements will be sent to the
 61:         appropriate processor during matrix assembly).
 62:       - Always specify global rows and columns of matrix entries.

 64:      Note: this uses the less common natural ordering that orders first
 65:      all the unknowns for x = h then for x = 2h etc; Hence you see J = Ii +- n
 66:      instead of J = I +- m as you might expect. The more standard ordering
 67:      would first do all variables for y = h, then y = 2h etc.

 69:    */
 70:   for (Ii=Istart; Ii<Iend; Ii++) {
 71:     v = -1.0; i = Ii/n; j = Ii - i*n;
 72:     if (i>0)   {J = Ii - n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
 73:     if (i<m-1) {J = Ii + n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
 74:     if (j>0)   {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
 75:     if (j<n-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
 76:     v = 4.0; MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);
 77:   }

 79:   /*
 80:      Assemble matrix, using the 2-step process:
 81:        MatAssemblyBegin(), MatAssemblyEnd()
 82:      Computations can be done while messages are in transition
 83:      by placing code between these two statements.
 84:   */
 85:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 86:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 88:   /* A is symmetric. Set symmetric flag to enable ICC/Cholesky preconditioner */
 89:   MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);

 91:   /*
 92:      Create parallel vectors.
 93:       - We form 1 vector from scratch and then duplicate as needed.
 94:       - When using VecCreate(), VecSetSizes and VecSetFromOptions()
 95:         in this example, we specify only the
 96:         vector's global dimension; the parallel partitioning is determined
 97:         at runtime.
 98:       - When solving a linear system, the vectors and matrices MUST
 99:         be partitioned accordingly.  PETSc automatically generates
100:         appropriately partitioned matrices and vectors when MatCreate()
101:         and VecCreate() are used with the same communicator.
102:       - The user can alternatively specify the local vector and matrix
103:         dimensions when more sophisticated partitioning is needed
104:         (replacing the PETSC_DECIDE argument in the VecSetSizes() statement
105:         below).
106:   */
107:   VecCreate(PETSC_COMM_WORLD,&u);
108:   VecSetSizes(u,PETSC_DECIDE,m*n);
109:   VecSetFromOptions(u);
110:   VecDuplicate(u,&b);
111:   VecDuplicate(b,&x);

113:   /*
114:      Set exact solution; then compute right-hand-side vector.
115:      By default we use an exact solution of a vector with all
116:      elements of 1.0;
117:   */
118:   VecSet(u,1.0);
119:   MatMult(A,u,b);

121:   /*
122:      View the exact solution vector if desired
123:   */
124:   flg  = PETSC_FALSE;
125:   PetscOptionsGetBool(NULL,NULL,"-view_exact_sol",&flg,NULL);
126:   if (flg) VecView(u,PETSC_VIEWER_STDOUT_WORLD);

128:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129:                 Create the linear solver and set various options
130:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
131:   KSPCreate(PETSC_COMM_WORLD,&ksp);

133:   /*
134:      Set operators. Here the matrix that defines the linear system
135:      also serves as the preconditioning matrix.
136:   */
137:   KSPSetOperators(ksp,A,A);

139:   /*
140:      Set linear solver defaults for this problem (optional).
141:      - By extracting the KSP and PC contexts from the KSP context,
142:        we can then directly call any KSP and PC routines to set
143:        various options.
144:      - The following two statements are optional; all of these
145:        parameters could alternatively be specified at runtime via
146:        KSPSetFromOptions().  All of these defaults can be
147:        overridden at runtime, as indicated below.
148:   */
149:   KSPSetTolerances(ksp,1.e-2/((m+1)*(n+1)),1.e-50,PETSC_DEFAULT,PETSC_DEFAULT);

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

160:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
161:                       Solve the linear system
162:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

164:   KSPSolve(ksp,b,x);

166:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167:                       Check the solution and clean up
168:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
169:   VecAXPY(x,-1.0,u);
170:   VecNorm(x,NORM_2,&norm);
171:   KSPGetIterationNumber(ksp,&its);

173:   /*
174:      Print convergence information.  PetscPrintf() produces a single
175:      print statement from all processes that share a communicator.
176:      An alternative is PetscFPrintf(), which prints to a file.
177:   */
178:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm,its);

180:   /*
181:      Free work space.  All PETSc objects should be destroyed when they
182:      are no longer needed.
183:   */
184:   KSPDestroy(&ksp);
185:   VecDestroy(&u));  PetscCall(VecDestroy(&x);
186:   VecDestroy(&b));  PetscCall(MatDestroy(&A);

188:   /*
189:      Always call PetscFinalize() before exiting a program.  This routine
190:        - finalizes the PETSc libraries as well as MPI
191:        - provides summary and diagnostic information if certain runtime
192:          options are chosen (e.g., -log_view).
193:   */
194:   PetscFinalize();
195:   return 0;
196: }

198: /*TEST

200:    build:
201:       requires: !single

203:    test:
204:       suffix: chebyest_1
205:       args: -m 80 -n 80 -ksp_pc_side right -pc_type ksp -ksp_ksp_type chebyshev -ksp_ksp_max_it 5 -ksp_ksp_chebyshev_esteig 0.9,0,0,1.1 -ksp_monitor_short

207:    test:
208:       suffix: chebyest_2
209:       args: -m 80 -n 80 -ksp_pc_side right -pc_type ksp -ksp_ksp_type chebyshev -ksp_ksp_max_it 5 -ksp_ksp_chebyshev_esteig 0.9,0,0,1.1 -ksp_esteig_ksp_type cg -ksp_monitor_short

211:    test:
212:       args: -ksp_monitor_short -m 5 -n 5 -ksp_gmres_cgs_refinement_type refine_always

214:    test:
215:       suffix: 2
216:       nsize: 2
217:       args: -ksp_monitor_short -m 5 -n 5 -ksp_gmres_cgs_refinement_type refine_always

219:    test:
220:       suffix: 3
221:       args: -pc_type sor -pc_sor_symmetric -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always

223:    test:
224:       suffix: 4
225:       args: -pc_type eisenstat -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always

227:    test:
228:       suffix: 5
229:       nsize: 2
230:       args: -ksp_monitor_short -m 5 -n 5 -mat_view draw -ksp_gmres_cgs_refinement_type refine_always -nox
231:       output_file: output/ex2_2.out

233:    test:
234:       suffix: bjacobi
235:       nsize: 4
236:       args: -pc_type bjacobi -pc_bjacobi_blocks 1 -ksp_monitor_short -sub_pc_type jacobi -sub_ksp_type gmres

238:    test:
239:       suffix: bjacobi_2
240:       nsize: 4
241:       args: -pc_type bjacobi -pc_bjacobi_blocks 2 -ksp_monitor_short -sub_pc_type jacobi -sub_ksp_type gmres -ksp_view

243:    test:
244:       suffix: bjacobi_3
245:       nsize: 4
246:       args: -pc_type bjacobi -pc_bjacobi_blocks 4 -ksp_monitor_short -sub_pc_type jacobi -sub_ksp_type gmres

248:    test:
249:       suffix: qmrcgs
250:       args: -ksp_type qmrcgs -pc_type ilu
251:       output_file: output/ex2_fbcgs.out

253:    test:
254:       suffix: qmrcgs_2
255:       nsize: 3
256:       args: -ksp_type qmrcgs -pc_type bjacobi
257:       output_file: output/ex2_fbcgs_2.out

259:    test:
260:       suffix: fbcgs
261:       args: -ksp_type fbcgs -pc_type ilu

263:    test:
264:       suffix: fbcgs_2
265:       nsize: 3
266:       args: -ksp_type fbcgsr -pc_type bjacobi

268:    test:
269:       suffix: groppcg
270:       args: -ksp_monitor_short -ksp_type groppcg -m 9 -n 9

272:    test:
273:       suffix: mkl_pardiso_cholesky
274:       requires: mkl_pardiso
275:       args: -ksp_type preonly -pc_type cholesky -mat_type sbaij -pc_factor_mat_solver_type mkl_pardiso

277:    test:
278:       suffix: mkl_pardiso_lu
279:       requires: mkl_pardiso
280:       args: -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mkl_pardiso

282:    test:
283:       suffix: pipebcgs
284:       args: -ksp_monitor_short -ksp_type pipebcgs -m 9 -n 9

286:    test:
287:       suffix: pipecg
288:       args: -ksp_monitor_short -ksp_type pipecg -m 9 -n 9

290:    test:
291:       suffix: pipecgrr
292:       args: -ksp_monitor_short -ksp_type pipecgrr -m 9 -n 9

294:    test:
295:       suffix: pipecr
296:       args: -ksp_monitor_short -ksp_type pipecr -m 9 -n 9

298:    test:
299:       suffix: pipelcg
300:       args: -ksp_monitor_short -ksp_type pipelcg -m 9 -n 9 -pc_type none -ksp_pipelcg_pipel 2 -ksp_pipelcg_lmax 2
301:       filter: grep -v "sqrt breakdown in iteration"

303:    test:
304:       suffix: sell
305:       args: -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always -m 9 -n 9 -mat_type sell

307:    test:
308:       requires: mumps
309:       suffix: sell_mumps
310:       args: -ksp_type preonly -m 9 -n 12 -mat_type sell -pc_type lu -pc_factor_mat_solver_type mumps -pc_factor_mat_ordering_type natural

312:    test:
313:       suffix: telescope
314:       nsize: 4
315:       args: -m 100 -n 100 -ksp_converged_reason -pc_type telescope -pc_telescope_reduction_factor 4 -telescope_pc_type bjacobi

317:    test:
318:       suffix: umfpack
319:       requires: suitesparse
320:       args: -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type umfpack

322:    test:
323:       suffix: spqr
324:       requires: suitesparse
325:       args: -ksp_type preonly -pc_type qr -pc_factor_mat_solver_type spqr

327:    test:
328:      suffix: pc_symmetric
329:      args: -m 10 -n 9 -ksp_converged_reason -ksp_type gmres -ksp_pc_side symmetric -pc_type cholesky

331:    test:
332:       suffix: pipeprcg
333:       args: -ksp_monitor_short -ksp_type pipeprcg -m 9 -n 9

335:    test:
336:       suffix: pipeprcg_rcw
337:       args: -ksp_monitor_short -ksp_type pipeprcg -recompute_w false -m 9 -n 9

339:    test:
340:       suffix: pipecg2
341:       args: -ksp_monitor_short -ksp_type pipecg2 -m 9 -n 9 -ksp_norm_type {{preconditioned unpreconditioned natural}}

343:    test:
344:       suffix: pipecg2_2
345:       nsize: 4
346:       args: -ksp_monitor_short -ksp_type pipecg2 -m 15 -n 9 -ksp_norm_type {{preconditioned unpreconditioned natural}}

348:    test:
349:       suffix: hpddm
350:       nsize: 4
351:       requires: hpddm
352:       filter: sed -e "s/ iterations 9/ iterations 8/g"
353:       args: -ksp_converged_reason -ksp_type hpddm -ksp_hpddm_precision {{single double}shared output} -ksp_pc_side {{left right}shared output}

355:  TEST*/