Actual source code: ex52.c
petsc-3.5.4 2015-05-23
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 and SUPERLU \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>
14: int main(int argc,char **args)
15: {
16: Vec x,b,u; /* approx solution, RHS, exact solution */
17: Mat A; /* linear system matrix */
18: KSP ksp; /* linear solver context */
19: PetscRandom rctx; /* random number generator context */
20: PetscReal norm; /* norm of solution error */
21: PetscInt i,j,Ii,J,Istart,Iend,m = 8,n = 7,its;
23: PetscBool flg,flg_ilu,flg_ch;
24: PetscScalar v;
25: PetscMPIInt rank;
26: #if defined(PETSC_USE_LOG)
27: PetscLogStage stage;
28: #endif
30: PetscInitialize(&argc,&args,(char*)0,help);
31: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
32: PetscOptionsGetInt(NULL,"-m",&m,NULL);
33: PetscOptionsGetInt(NULL,"-n",&n,NULL);
34: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
35: Compute the matrix and right-hand-side vector that define
36: the linear system, Ax = b.
37: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
38: MatCreate(PETSC_COMM_WORLD,&A);
39: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
40: MatSetFromOptions(A);
41: MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);
42: MatSeqAIJSetPreallocation(A,5,NULL);
43: MatSetUp(A);
45: /*
46: Currently, all PETSc parallel matrix formats are partitioned by
47: contiguous chunks of rows across the processors. Determine which
48: rows of the matrix are locally owned.
49: */
50: MatGetOwnershipRange(A,&Istart,&Iend);
52: /*
53: Set matrix elements for the 2-D, five-point stencil in parallel.
54: - Each processor needs to insert only elements that it owns
55: locally (but any non-local elements will be sent to the
56: appropriate processor during matrix assembly).
57: - Always specify global rows and columns of matrix entries.
59: Note: this uses the less common natural ordering that orders first
60: all the unknowns for x = h then for x = 2h etc; Hence you see J = Ii +- n
61: instead of J = I +- m as you might expect. The more standard ordering
62: would first do all variables for y = h, then y = 2h etc.
64: */
65: PetscLogStageRegister("Assembly", &stage);
66: PetscLogStagePush(stage);
67: for (Ii=Istart; Ii<Iend; Ii++) {
68: v = -1.0; i = Ii/n; j = Ii - i*n;
69: if (i>0) {J = Ii - n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
70: if (i<m-1) {J = Ii + n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
71: if (j>0) {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
72: if (j<n-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
73: v = 4.0; MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES);
74: }
76: /*
77: Assemble matrix, using the 2-step process:
78: MatAssemblyBegin(), MatAssemblyEnd()
79: Computations can be done while messages are in transition
80: by placing code between these two statements.
81: */
82: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
83: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
84: PetscLogStagePop();
86: /* A is symmetric. Set symmetric flag to enable ICC/Cholesky preconditioner */
87: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
89: /*
90: Create parallel vectors.
91: - We form 1 vector from scratch and then duplicate as needed.
92: - When using VecCreate(), VecSetSizes and VecSetFromOptions()
93: in this example, we specify only the
94: vector's global dimension; the parallel partitioning is determined
95: at runtime.
96: - When solving a linear system, the vectors and matrices MUST
97: be partitioned accordingly. PETSc automatically generates
98: appropriately partitioned matrices and vectors when MatCreate()
99: and VecCreate() are used with the same communicator.
100: - The user can alternatively specify the local vector and matrix
101: dimensions when more sophisticated partitioning is needed
102: (replacing the PETSC_DECIDE argument in the VecSetSizes() statement
103: below).
104: */
105: VecCreate(PETSC_COMM_WORLD,&u);
106: VecSetSizes(u,PETSC_DECIDE,m*n);
107: VecSetFromOptions(u);
108: VecDuplicate(u,&b);
109: VecDuplicate(b,&x);
111: /*
112: Set exact solution; then compute right-hand-side vector.
113: By default we use an exact solution of a vector with all
114: elements of 1.0; Alternatively, using the runtime option
115: -random_sol forms a solution vector with random components.
116: */
117: PetscOptionsGetBool(NULL,"-random_exact_sol",&flg,NULL);
118: if (flg) {
119: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
120: PetscRandomSetFromOptions(rctx);
121: VecSetRandom(u,rctx);
122: PetscRandomDestroy(&rctx);
123: } else {
124: VecSet(u,1.0);
125: }
126: MatMult(A,u,b);
128: /*
129: View the exact solution vector if desired
130: */
131: flg = PETSC_FALSE;
132: PetscOptionsGetBool(NULL,"-view_exact_sol",&flg,NULL);
133: if (flg) {VecView(u,PETSC_VIEWER_STDOUT_WORLD);}
135: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136: Create the linear solver and set various options
137: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
139: /*
140: Create linear solver context
141: */
142: KSPCreate(PETSC_COMM_WORLD,&ksp);
143: KSPSetOperators(ksp,A,A);
145: /*
146: Example of how to use external package MUMPS
147: Note: runtime options
148: '-ksp_type preonly -pc_type lu -pc_factor_mat_solver_package mumps -mat_mumps_icntl_7 2 -mat_mumps_icntl_1 0.0'
149: are equivalent to these procedural calls
150: */
151: #if defined(PETSC_HAVE_MUMPS)
152: Mat F;
153: flg = PETSC_FALSE;
154: flg_ch = PETSC_FALSE;
155: PetscOptionsGetBool(NULL,"-use_mumps_lu",&flg,NULL);
156: PetscOptionsGetBool(NULL,"-use_mumps_ch",&flg_ch,NULL);
157: if (flg || flg_ch) {
158: KSPSetType(ksp,KSPPREONLY);
159: PC pc;
160: PetscInt ival,icntl;
161: PetscReal val;
162: KSPGetPC(ksp,&pc);
163: if (flg) {
164: PCSetType(pc,PCLU);
165: } else if (flg_ch) {
166: MatSetOption(A,MAT_SPD,PETSC_TRUE); /* set MUMPS id%SYM=1 */
167: PCSetType(pc,PCCHOLESKY);
168: }
169: PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS);
170: PCFactorSetUpMatSolverPackage(pc); /* call MatGetFactor() to create F */
171: PCFactorGetMatrix(pc,&F);
173: /* sequential ordering */
174: icntl = 7; ival = 2;
175: MatMumpsSetIcntl(F,icntl,ival);
177: /* threshhold for row pivot detection */
178: MatMumpsSetIcntl(F,24,1);
179: icntl = 3; val = 1.e-6;
180: MatMumpsSetCntl(F,icntl,val);
182: /* compute determinant of A */
183: MatMumpsSetIcntl(F,33,1);
184: }
185: #endif
187: /*
188: Example of how to use external package SuperLU
189: Note: runtime options
190: '-ksp_type preonly -pc_type ilu -pc_factor_mat_solver_package superlu -mat_superlu_ilu_droptol 1.e-8'
191: are equivalent to these procedual calls
192: */
193: #if defined(PETSC_HAVE_SUPERLU)
194: flg_ilu = PETSC_FALSE;
195: flg = PETSC_FALSE;
196: PetscOptionsGetBool(NULL,"-use_superlu_lu",&flg,NULL);
197: PetscOptionsGetBool(NULL,"-use_superlu_ilu",&flg_ilu,NULL);
198: if (flg || flg_ilu) {
199: KSPSetType(ksp,KSPPREONLY);
200: PC pc;
201: Mat F;
202: KSPGetPC(ksp,&pc);
203: if (flg) {
204: PCSetType(pc,PCLU);
205: } else if (flg_ilu) {
206: PCSetType(pc,PCILU);
207: }
208: PCFactorSetMatSolverPackage(pc,MATSOLVERSUPERLU);
209: PCFactorSetUpMatSolverPackage(pc); /* call MatGetFactor() to create F */
210: PCFactorGetMatrix(pc,&F);
211: MatSuperluSetILUDropTol(F,1.e-8);
212: }
213: #endif
215: /*
216: Example of how to use procedural calls that are equivalent to
217: '-ksp_type preonly -pc_type lu/ilu -pc_factor_mat_solver_package petsc'
218: */
219: flg = PETSC_FALSE;
220: flg_ilu = PETSC_FALSE;
221: flg_ch = PETSC_FALSE;
222: PetscOptionsGetBool(NULL,"-use_petsc_lu",&flg,NULL);
223: PetscOptionsGetBool(NULL,"-use_petsc_ilu",&flg_ilu,NULL);
224: PetscOptionsGetBool(NULL,"-use_petsc_ch",&flg_ch,NULL);
225: if (flg || flg_ilu || flg_ch) {
226: KSPSetType(ksp,KSPPREONLY);
227: PC pc;
228: Mat F;
229: KSPGetPC(ksp,&pc);
230: if (flg) {
231: PCSetType(pc,PCLU);
232: } else if (flg_ilu) {
233: PCSetType(pc,PCILU);
234: } else if (flg_ch) {
235: PCSetType(pc,PCCHOLESKY);
236: }
237: PCFactorSetMatSolverPackage(pc,MATSOLVERPETSC);
238: PCFactorSetUpMatSolverPackage(pc); /* call MatGetFactor() to create F */
239: PCFactorGetMatrix(pc,&F);
241: /* Test MatGetDiagonal() */
242: Vec diag;
243: KSPSetUp(ksp);
244: MatView(F,PETSC_VIEWER_STDOUT_WORLD);
245: VecDuplicate(x,&diag);
246: MatGetDiagonal(F,diag);
247: VecView(diag,PETSC_VIEWER_STDOUT_WORLD);
248: VecDestroy(&diag);
249: }
251: KSPSetFromOptions(ksp);
253: /* Get info from matrix factors */
254: KSPSetUp(ksp);
256: #if defined(PETSC_HAVE_MUMPS)
257: if (flg || flg_ch) {
258: PetscInt icntl,infog34;
259: PetscReal cntl,rinfo12,rinfo13;
260: icntl = 3;
261: MatMumpsGetCntl(F,icntl,&cntl);
262:
263: /* compute determinant */
264: if (!rank) {
265: MatMumpsGetInfog(F,34,&infog34);
266: MatMumpsGetRinfog(F,12,&rinfo12);
267: MatMumpsGetRinfog(F,13,&rinfo13);
268: PetscPrintf(PETSC_COMM_SELF," Mumps row pivot threshhold = %g\n",cntl);
269: PetscPrintf(PETSC_COMM_SELF," Mumps determinant = (%g, %g) * 2^%D \n",rinfo12,rinfo13,infog34);
270: }
271: }
272: #endif
274: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
275: Solve the linear system
276: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
277: KSPSolve(ksp,b,x);
279: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
280: Check solution and clean up
281: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
282: VecAXPY(x,-1.0,u);
283: VecNorm(x,NORM_2,&norm);
284: KSPGetIterationNumber(ksp,&its);
286: /*
287: Print convergence information. PetscPrintf() produces a single
288: print statement from all processes that share a communicator.
289: An alternative is PetscFPrintf(), which prints to a file.
290: */
291: if (norm < 1.e-12) {
292: PetscPrintf(PETSC_COMM_WORLD,"Norm of error < 1.e-12 iterations %D\n",norm,its);
293: } else {
294: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm,its);
295: }
297: /*
298: Free work space. All PETSc objects should be destroyed when they
299: are no longer needed.
300: */
301: KSPDestroy(&ksp);
302: VecDestroy(&u); VecDestroy(&x);
303: VecDestroy(&b); MatDestroy(&A);
305: /*
306: Always call PetscFinalize() before exiting a program. This routine
307: - finalizes the PETSc libraries as well as MPI
308: - provides summary and diagnostic information if certain runtime
309: options are chosen (e.g., -log_summary).
310: */
311: PetscFinalize();
312: return 0;
313: }