Actual source code: ex52.c
petsc-3.9.4 2018-09-11
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*/