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