Actual source code: ex52.c

petsc-master 2018-05-23
Report Typos and Errors

  2: static char help[] = "Solves a linear system in parallel with KSP. Modified from ex2.c \n\
  3:                       Illustrate how to use external packages MUMPS, SUPERLU and STRUMPACK \n\
  4: Input parameters include:\n\
  5:   -random_exact_sol : use a random exact solution vector\n\
  6:   -view_exact_sol   : write exact solution vector to stdout\n\
  7:   -m <mesh_x>       : number of mesh points in x-direction\n\
  8:   -n <mesh_n>       : number of mesh points in y-direction\n\n";

 10:  #include <petscksp.h>

 12: int main(int argc,char **args)
 13: {
 14:   Vec            x,b,u;    /* approx solution, RHS, exact solution */
 15:   Mat            A,F;
 16:   KSP            ksp;      /* linear solver context */
 17:   PC             pc;
 18:   PetscRandom    rctx;     /* random number generator context */
 19:   PetscReal      norm;     /* norm of solution error */
 20:   PetscInt       i,j,Ii,J,Istart,Iend,m = 8,n = 7,its;
 22:   PetscBool      flg=PETSC_FALSE,flg_ilu=PETSC_FALSE,flg_ch=PETSC_FALSE;
 23: #if defined(PETSC_HAVE_MUMPS)
 24:   PetscBool      flg_mumps=PETSC_FALSE,flg_mumps_ch=PETSC_FALSE;
 25: #endif
 26: #if defined(PETSC_HAVE_SUPERLU) || defined(PETSC_HAVE_SUPERLU_DIST)
 27:   PetscBool      flg_superlu=PETSC_FALSE;
 28: #endif
 29: #if defined(PETSC_HAVE_STRUMPACK)
 30:   PetscBool      flg_strumpack=PETSC_FALSE;
 31: #endif
 32:   PetscScalar    v;
 33:   PetscMPIInt    rank,size;
 34: #if defined(PETSC_USE_LOG)
 35:   PetscLogStage  stage;
 36: #endif

 38:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 39:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 40:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 41:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
 42:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 43:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 44:          Compute the matrix and right-hand-side vector that define
 45:          the linear system, Ax = b.
 46:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 47:   MatCreate(PETSC_COMM_WORLD,&A);
 48:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
 49:   MatSetFromOptions(A);
 50:   MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);
 51:   MatSeqAIJSetPreallocation(A,5,NULL);
 52:   MatSetUp(A);

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

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

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

 73:    */
 74:   PetscLogStageRegister("Assembly", &stage);
 75:   PetscLogStagePush(stage);
 76:   for (Ii=Istart; Ii<Iend; Ii++) {
 77:     v = -1.0; i = Ii/n; j = Ii - i*n;
 78:     if (i>0)   {J = Ii - n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
 79:     if (i<m-1) {J = Ii + n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
 80:     if (j>0)   {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
 81:     if (j<n-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
 82:     v = 4.0; MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES);
 83:   }

 85:   /*
 86:      Assemble matrix, using the 2-step process:
 87:        MatAssemblyBegin(), MatAssemblyEnd()
 88:      Computations can be done while messages are in transition
 89:      by placing code between these two statements.
 90:   */
 91:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 92:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 93:   PetscLogStagePop();

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

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

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

137:   /*
138:      View the exact solution vector if desired
139:   */
140:   flg  = PETSC_FALSE;
141:   PetscOptionsGetBool(NULL,NULL,"-view_exact_sol",&flg,NULL);
142:   if (flg) {VecView(u,PETSC_VIEWER_STDOUT_WORLD);}

144:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145:                 Create the linear solver and set various options
146:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

148:   /*
149:      Create linear solver context
150:   */
151:   KSPCreate(PETSC_COMM_WORLD,&ksp);
152:   KSPSetOperators(ksp,A,A);

154:   /*
155:     Example of how to use external package MUMPS
156:     Note: runtime options
157:           '-ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps -mat_mumps_icntl_7 2 -mat_mumps_icntl_1 0.0'
158:           are equivalent to these procedural calls
159:   */
160: #if defined(PETSC_HAVE_MUMPS)
161:   flg_mumps    = PETSC_FALSE;
162:   flg_mumps_ch = PETSC_FALSE;
163:   PetscOptionsGetBool(NULL,NULL,"-use_mumps_lu",&flg_mumps,NULL);
164:   PetscOptionsGetBool(NULL,NULL,"-use_mumps_ch",&flg_mumps_ch,NULL);
165:   if (flg_mumps || flg_mumps_ch) {
166:     KSPSetType(ksp,KSPPREONLY);
167:     PetscInt  ival,icntl;
168:     PetscReal val;
169:     KSPGetPC(ksp,&pc);
170:     if (flg_mumps) {
171:       PCSetType(pc,PCLU);
172:     } else if (flg_mumps_ch) {
173:       MatSetOption(A,MAT_SPD,PETSC_TRUE); /* set MUMPS id%SYM=1 */
174:       PCSetType(pc,PCCHOLESKY);
175:     }
176:     PCFactorSetMatSolverType(pc,MATSOLVERMUMPS);
177:     PCFactorSetUpMatSolverType(pc); /* call MatGetFactor() to create F */
178:     PCFactorGetMatrix(pc,&F);

180:     /* sequential ordering */
181:     icntl = 7; ival = 2;
182:     MatMumpsSetIcntl(F,icntl,ival);

184:     /* threshhold for row pivot detection */
185:     MatMumpsSetIcntl(F,24,1);
186:     icntl = 3; val = 1.e-6;
187:     MatMumpsSetCntl(F,icntl,val);

189:     /* compute determinant of A */
190:     MatMumpsSetIcntl(F,33,1);
191:   }
192: #endif

194:   /*
195:     Example of how to use external package SuperLU
196:     Note: runtime options
197:           '-ksp_type preonly -pc_type ilu -pc_factor_mat_solver_type superlu -mat_superlu_ilu_droptol 1.e-8'
198:           are equivalent to these procedual calls
199:   */
200: #if defined(PETSC_HAVE_SUPERLU) || defined(PETSC_HAVE_SUPERLU_DIST)
201:   flg_ilu     = PETSC_FALSE;
202:   flg_superlu = PETSC_FALSE;
203:   PetscOptionsGetBool(NULL,NULL,"-use_superlu_lu",&flg_superlu,NULL);
204:   PetscOptionsGetBool(NULL,NULL,"-use_superlu_ilu",&flg_ilu,NULL);
205:   if (flg_superlu || flg_ilu) {
206:     KSPSetType(ksp,KSPPREONLY);
207:     KSPGetPC(ksp,&pc);
208:     if (flg_superlu) {
209:       PCSetType(pc,PCLU);
210:     } else if (flg_ilu) {
211:       PCSetType(pc,PCILU);
212:     }
213:     if (size == 1) {
214: #if !defined(PETSC_HAVE_SUPERLU)
215:       SETERRQ(PETSC_COMM_WORLD,1,"This test requires SUPERLU");
216: #else
217:       PCFactorSetMatSolverType(pc,MATSOLVERSUPERLU);
218: #endif
219:     } else {
220: #if !defined(PETSC_HAVE_SUPERLU_DIST)
221:       SETERRQ(PETSC_COMM_WORLD,1,"This test requires SUPERLU_DIST");
222: #else
223:       PCFactorSetMatSolverType(pc,MATSOLVERSUPERLU_DIST);
224: #endif
225:     }
226:     PCFactorSetUpMatSolverType(pc); /* call MatGetFactor() to create F */
227:     PCFactorGetMatrix(pc,&F);
228: #if defined(PETSC_HAVE_SUPERLU)
229:     if (size == 1) {
230:       MatSuperluSetILUDropTol(F,1.e-8);
231:     }
232: #endif
233:   }
234: #endif


237:   /*
238:     Example of how to use external package STRUMPACK
239:     Note: runtime options
240:           '-pc_type lu/ilu \
241:            -pc_factor_mat_solver_type strumpack \
242:            -mat_strumpack_reordering METIS \
243:            -mat_strumpack_colperm 0 \
244:            -mat_strumpack_hss_rel_tol 1.e-3 \
245:            -mat_strumpack_hss_min_sep_size 50 \
246:            -mat_strumpack_max_rank 100 \
247:            -mat_strumpack_leaf_size 4'
248:        are equivalent to these procedural calls

250:     We refer to the STRUMPACK-sparse manual, section 5, for more info on
251:     how to tune the preconditioner.
252:   */
253: #if defined(PETSC_HAVE_STRUMPACK)
254:   flg_ilu       = PETSC_FALSE;
255:   flg_strumpack = PETSC_FALSE;
256:   PetscOptionsGetBool(NULL,NULL,"-use_strumpack_lu",&flg_strumpack,NULL);
257:   PetscOptionsGetBool(NULL,NULL,"-use_strumpack_ilu",&flg_ilu,NULL);
258:   if (flg_strumpack || flg_ilu) {
259:     KSPSetType(ksp,KSPPREONLY);
260:     KSPGetPC(ksp,&pc);
261:     if (flg_strumpack) {
262:       PCSetType(pc,PCLU);
263:     } else if (flg_ilu) {
264:       PCSetType(pc,PCILU);
265:     }
266: #if !defined(PETSC_HAVE_STRUMPACK)
267:     SETERRQ(PETSC_COMM_WORLD,1,"This test requires STRUMPACK");
268: #endif
269:     PCFactorSetMatSolverType(pc,MATSOLVERSTRUMPACK);
270:     PCFactorSetUpMatSolverType(pc); /* call MatGetFactor() to create F */
271:     PCFactorGetMatrix(pc,&F);
272: #if defined(PETSC_HAVE_STRUMPACK)
273:     /* Set the fill-reducing reordering.                              */
274:     MatSTRUMPACKSetReordering(F,MAT_STRUMPACK_METIS);
275:     /* Since this is a simple discretization, the diagonal is always  */
276:     /* nonzero, and there is no need for the extra MC64 permutation.  */
277:     MatSTRUMPACKSetColPerm(F,PETSC_FALSE);
278:     /* The compression tolerance used when doing low-rank compression */
279:     /* in the preconditioner. This is problem specific!               */
280:     MatSTRUMPACKSetHSSRelTol(F,1.e-3);
281:     /* Set minimum matrix size for HSS compression to 15 in order to  */
282:     /* demonstrate preconditioner on small problems. For performance  */
283:     /* a value of say 500 is better.                                  */
284:     MatSTRUMPACKSetHSSMinSepSize(F,15);
285:     /* You can further limit the fill in the preconditioner by        */
286:     /* setting a maximum rank                                         */
287:     MatSTRUMPACKSetHSSMaxRank(F,100);
288:     /* Set the size of the diagonal blocks (the leafs) in the HSS     */
289:     /* approximation. The default value should be better for real     */
290:     /* problems. This is mostly for illustration on a small problem.  */
291:     MatSTRUMPACKSetHSSLeafSize(F,4);
292: #endif
293:   }
294: #endif

296:   /*
297:     Example of how to use procedural calls that are equivalent to
298:           '-ksp_type preonly -pc_type lu/ilu -pc_factor_mat_solver_type petsc'
299:   */
300:   flg     = PETSC_FALSE;
301:   flg_ilu = PETSC_FALSE;
302:   flg_ch  = PETSC_FALSE;
303:   PetscOptionsGetBool(NULL,NULL,"-use_petsc_lu",&flg,NULL);
304:   PetscOptionsGetBool(NULL,NULL,"-use_petsc_ilu",&flg_ilu,NULL);
305:   PetscOptionsGetBool(NULL,NULL,"-use_petsc_ch",&flg_ch,NULL);
306:   if (flg || flg_ilu || flg_ch) {
307:     Vec diag;

309:     KSPSetType(ksp,KSPPREONLY);
310:     KSPGetPC(ksp,&pc);
311:     if (flg) {
312:       PCSetType(pc,PCLU);
313:     } else if (flg_ilu) {
314:       PCSetType(pc,PCILU);
315:     } else if (flg_ch) {
316:       PCSetType(pc,PCCHOLESKY);
317:     }
318:     PCFactorSetMatSolverType(pc,MATSOLVERPETSC);
319:     PCFactorSetUpMatSolverType(pc); /* call MatGetFactor() to create F */
320:     PCFactorGetMatrix(pc,&F);
321: 
322:     /* Test MatGetDiagonal() */
323:     KSPSetUp(ksp);
324:     VecDuplicate(x,&diag);
325:     MatGetDiagonal(F,diag);
326:     /* VecView(diag,PETSC_VIEWER_STDOUT_WORLD); */
327:     VecDestroy(&diag);
328:   }

330:   KSPSetFromOptions(ksp);

332:   /* Get info from matrix factors */
333:   KSPSetUp(ksp);

335: #if defined(PETSC_HAVE_MUMPS)
336:   if (flg_mumps || flg_mumps_ch) {
337:     PetscInt  icntl,infog34;
338:     PetscReal cntl,rinfo12,rinfo13;
339:     icntl = 3;
340:     MatMumpsGetCntl(F,icntl,&cntl);
341: 
342:     /* compute determinant */
343:     if (!rank) {
344:       MatMumpsGetInfog(F,34,&infog34);
345:       MatMumpsGetRinfog(F,12,&rinfo12);
346:       MatMumpsGetRinfog(F,13,&rinfo13);
347:       PetscPrintf(PETSC_COMM_SELF,"  Mumps row pivot threshhold = %g\n",cntl);
348:       PetscPrintf(PETSC_COMM_SELF,"  Mumps determinant = (%g, %g) * 2^%D \n",(double)rinfo12,(double)rinfo13,infog34);
349:     }
350:   }
351: #endif

353:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
354:                       Solve the linear system
355:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
356:   KSPSolve(ksp,b,x);

358:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
359:                       Check solution and clean up
360:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
361:   VecAXPY(x,-1.0,u);
362:   VecNorm(x,NORM_2,&norm);
363:   KSPGetIterationNumber(ksp,&its);

365:   /*
366:      Print convergence information.  PetscPrintf() produces a single
367:      print statement from all processes that share a communicator.
368:      An alternative is PetscFPrintf(), which prints to a file.
369:   */
370:   if (norm < 1.e-12) {
371:     PetscPrintf(PETSC_COMM_WORLD,"Norm of error < 1.e-12 iterations %D\n",its);
372:   } else {
373:     PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm,its);
374:  }

376:   /*
377:      Free work space.  All PETSc objects should be destroyed when they
378:      are no longer needed.
379:   */
380:   KSPDestroy(&ksp);
381:   VecDestroy(&u);  VecDestroy(&x);
382:   VecDestroy(&b);  MatDestroy(&A);

384:   /*
385:      Always call PetscFinalize() before exiting a program.  This routine
386:        - finalizes the PETSc libraries as well as MPI
387:        - provides summary and diagnostic information if certain runtime
388:          options are chosen (e.g., -log_view).
389:   */
390:   PetscFinalize();
391:   return ierr;
392: }


395: /*TEST

397:    test:
398:       requires: mumps
399:       args: -use_petsc_lu
400:       output_file: output/ex52_2.out

402:    test:
403:       suffix: mumps
404:       nsize: 3
405:       requires: mumps
406:       args: -use_mumps_lu
407:       output_file: output/ex52_1.out

409:    test:
410:       suffix: mumps_2
411:       nsize: 3
412:       requires: mumps
413:       args: -use_mumps_ch
414:       output_file: output/ex52_1.out

416:    test:
417:       suffix: mumps_3
418:       nsize: 3
419:       requires: mumps
420:       args: -use_mumps_ch -mat_type sbaij
421:       output_file: output/ex52_1.out

423:    test:
424:       suffix: strumpack
425:       requires: strumpack
426:       args: -use_strumpack_lu 
427:       output_file: output/ex52_3.out

429:    test:
430:       suffix: strumpack_2
431:       nsize: 2
432:       requires: strumpack
433:       args: -use_strumpack_lu
434:       output_file: output/ex52_3.out

436:    test:
437:       suffix: strumpack_ilu
438:       requires: strumpack
439:       args: -use_strumpack_ilu
440:       output_file: output/ex52_3.out

442:    test:
443:       suffix: strumpack_ilu_2
444:       nsize: 2
445:       requires: strumpack
446:       args: -use_strumpack_ilu
447:       output_file: output/ex52_3.out

449:    test:
450:       suffix: superlu
451:       requires: superlu superlu_dist
452:       args: -use_superlu_lu
453:       output_file: output/ex52_2.out

455:    test:
456:       suffix: superlu_dist
457:       nsize: 2
458:       requires: superlu superlu_dist
459:       args: -use_superlu_lu
460:       output_file: output/ex52_2.out

462:    test:
463:       suffix: superlu_ilu
464:       requires: superlu superlu_dist
465:       args: -use_superlu_ilu
466:       output_file: output/ex52_2.out

468: TEST*/