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