Actual source code: ex52.c

petsc-master 2017-12-15
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_package 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:     PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS);
177:     PCFactorSetUpMatSolverPackage(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_package 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: #endif
217:       PCFactorSetMatSolverPackage(pc,MATSOLVERSUPERLU);
218:     } else {
219: #if !defined(PETSC_HAVE_SUPERLU_DIST)
220:       SETERRQ(PETSC_COMM_WORLD,1,"This test requires SUPERLU_DIST");
221: #endif
222:       PCFactorSetMatSolverPackage(pc,MATSOLVERSUPERLU_DIST);
223:     }
224:     PCFactorSetUpMatSolverPackage(pc); /* call MatGetFactor() to create F */
225:     PCFactorGetMatrix(pc,&F);
226: #if defined(PETSC_HAVE_SUPERLU)
227:     if (size == 1) {
228:       MatSuperluSetILUDropTol(F,1.e-8);
229:     }
230: #endif
231:   }
232: #endif


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

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

294:   /*
295:     Example of how to use procedural calls that are equivalent to
296:           '-ksp_type preonly -pc_type lu/ilu -pc_factor_mat_solver_package petsc'
297:   */
298:   flg     = PETSC_FALSE;
299:   flg_ilu = PETSC_FALSE;
300:   flg_ch  = PETSC_FALSE;
301:   PetscOptionsGetBool(NULL,NULL,"-use_petsc_lu",&flg,NULL);
302:   PetscOptionsGetBool(NULL,NULL,"-use_petsc_ilu",&flg_ilu,NULL);
303:   PetscOptionsGetBool(NULL,NULL,"-use_petsc_ch",&flg_ch,NULL);
304:   if (flg || flg_ilu || flg_ch) {
305:     KSPSetType(ksp,KSPPREONLY);
306:     KSPGetPC(ksp,&pc);
307:     if (flg) {
308:       PCSetType(pc,PCLU);
309:     } else if (flg_ilu) {
310:       PCSetType(pc,PCILU);
311:     } else if (flg_ch) {
312:       PCSetType(pc,PCCHOLESKY);
313:     }
314:     PCFactorSetMatSolverPackage(pc,MATSOLVERPETSC);
315:     PCFactorSetUpMatSolverPackage(pc); /* call MatGetFactor() to create F */
316:     PCFactorGetMatrix(pc,&F);
317: 
318:     /* Test MatGetDiagonal() */
319:     Vec diag;
320:     KSPSetUp(ksp);
321:     VecDuplicate(x,&diag);
322:     MatGetDiagonal(F,diag);
323:     /* VecView(diag,PETSC_VIEWER_STDOUT_WORLD); */
324:     VecDestroy(&diag);
325:   }

327:   KSPSetFromOptions(ksp);

329:   /* Get info from matrix factors */
330:   KSPSetUp(ksp);

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

350:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
351:                       Solve the linear system
352:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
353:   KSPSolve(ksp,b,x);

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

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

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

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