Actual source code: ex10.c
petsc-3.7.5 2017-01-01
2: static char help[] = "Reads a PETSc matrix and vector from a file and solves a linear system.\n\
3: This version first preloads and solves a small system, then loads \n\
4: another (larger) system and solves it as well. This example illustrates\n\
5: preloading of instructions with the smaller system so that more accurate\n\
6: performance monitoring can be done with the larger one (that actually\n\
7: is the system of interest). See the 'Performance Hints' chapter of the\n\
8: users manual for a discussion of preloading. Input parameters include\n\
9: -f0 <input_file> : first file to load (small system)\n\
10: -f1 <input_file> : second file to load (larger system)\n\n\
11: -nearnulldim <0> : number of vectors in the near-null space immediately following matrix\n\n\
12: -trans : solve transpose system instead\n\n";
13: /*
14: This code can be used to test PETSc interface to other packages.\n\
15: Examples of command line options: \n\
16: ./ex10 -f0 <datafile> -ksp_type preonly \n\
17: -help -ksp_view \n\
18: -num_numfac <num_numfac> -num_rhs <num_rhs> \n\
19: -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package superlu or superlu_dist or mumps \n\
20: -ksp_type preonly -pc_type cholesky -pc_factor_mat_solver_package mumps \n\
21: mpiexec -n <np> ./ex10 -f0 <datafile> -ksp_type cg -pc_type asm -pc_asm_type basic -sub_pc_type icc -mat_type sbaij
22: \n\n";
23: */
24: /*T
25: Concepts: KSP^solving a linear system
26: Processors: n
27: T*/
29: /*
30: Include "petscksp.h" so that we can use KSP solvers. Note that this file
31: automatically includes:
32: petscsys.h - base PETSc routines petscvec.h - vectors
33: petscmat.h - matrices
34: petscis.h - index sets petscksp.h - Krylov subspace methods
35: petscviewer.h - viewers petscpc.h - preconditioners
36: */
37: #include <petscksp.h>
41: int main(int argc,char **args)
42: {
43: KSP ksp; /* linear solver context */
44: Mat A; /* matrix */
45: Vec x,b,u; /* approx solution, RHS, exact solution */
46: PetscViewer fd; /* viewer */
47: char file[4][PETSC_MAX_PATH_LEN]; /* input file name */
48: PetscBool table =PETSC_FALSE,flg,trans=PETSC_FALSE,initialguess = PETSC_FALSE;
49: PetscBool outputSoln=PETSC_FALSE;
51: PetscInt its,num_numfac,m,n,M,nearnulldim = 0;
52: PetscReal norm;
53: PetscBool preload=PETSC_TRUE,isSymmetric,cknorm=PETSC_FALSE,initialguessfile = PETSC_FALSE;
54: PetscMPIInt rank;
55: char initialguessfilename[PETSC_MAX_PATH_LEN];
57: PetscInitialize(&argc,&args,(char*)0,help);
58: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
59: PetscOptionsGetBool(NULL,NULL,"-table",&table,NULL);
60: PetscOptionsGetBool(NULL,NULL,"-trans",&trans,NULL);
61: PetscOptionsGetBool(NULL,NULL,"-initialguess",&initialguess,NULL);
62: PetscOptionsGetBool(NULL,NULL,"-output_solution",&outputSoln,NULL);
63: PetscOptionsGetString(NULL,NULL,"-initialguessfilename",initialguessfilename,PETSC_MAX_PATH_LEN,&initialguessfile);
64: PetscOptionsGetInt(NULL,NULL,"-nearnulldim",&nearnulldim,NULL);
66: /*
67: Determine files from which we read the two linear systems
68: (matrix and right-hand-side vector).
69: */
70: PetscOptionsGetString(NULL,NULL,"-f",file[0],PETSC_MAX_PATH_LEN,&flg);
71: if (flg) {
72: PetscStrcpy(file[1],file[0]);
73: preload = PETSC_FALSE;
74: } else {
75: PetscOptionsGetString(NULL,NULL,"-f0",file[0],PETSC_MAX_PATH_LEN,&flg);
76: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f0 or -f option");
77: PetscOptionsGetString(NULL,NULL,"-f1",file[1],PETSC_MAX_PATH_LEN,&flg);
78: if (!flg) preload = PETSC_FALSE; /* don't bother with second system */
79: }
81: /* -----------------------------------------------------------
82: Beginning of linear solver loop
83: ----------------------------------------------------------- */
84: /*
85: Loop through the linear solve 2 times.
86: - The intention here is to preload and solve a small system;
87: then load another (larger) system and solve it as well.
88: This process preloads the instructions with the smaller
89: system so that more accurate performance monitoring (via
90: -log_summary) can be done with the larger one (that actually
91: is the system of interest).
92: */
93: PetscPreLoadBegin(preload,"Load system");
95: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
96: Load system
97: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99: /*
100: Open binary file. Note that we use FILE_MODE_READ to indicate
101: reading from this file.
102: */
103: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[PetscPreLoadIt],FILE_MODE_READ,&fd);
105: /*
106: Load the matrix and vector; then destroy the viewer.
107: */
108: MatCreate(PETSC_COMM_WORLD,&A);
109: MatSetFromOptions(A);
110: MatLoad(A,fd);
111: if (nearnulldim) {
112: MatNullSpace nullsp;
113: Vec *nullvecs;
114: PetscInt i;
115: PetscMalloc1(nearnulldim,&nullvecs);
116: for (i=0; i<nearnulldim; i++) {
117: VecCreate(PETSC_COMM_WORLD,&nullvecs[i]);
118: VecLoad(nullvecs[i],fd);
119: }
120: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_FALSE,nearnulldim,nullvecs,&nullsp);
121: MatSetNearNullSpace(A,nullsp);
122: for (i=0; i<nearnulldim; i++) {VecDestroy(&nullvecs[i]);}
123: PetscFree(nullvecs);
124: MatNullSpaceDestroy(&nullsp);
125: }
127: flg = PETSC_FALSE;
128: PetscOptionsGetString(NULL,NULL,"-rhs",file[2],PETSC_MAX_PATH_LEN,&flg);
129: VecCreate(PETSC_COMM_WORLD,&b);
130: if (flg) { /* rhs is stored in a separate file */
131: if (file[2][0] == '0' || file[2][0] == 0) {
132: PetscInt m;
133: PetscScalar one = 1.0;
134: PetscInfo(0,"Using vector of ones for RHS\n");
135: MatGetLocalSize(A,&m,NULL);
136: VecSetSizes(b,m,PETSC_DECIDE);
137: VecSetFromOptions(b);
138: VecSet(b,one);
139: } else {
140: PetscViewerDestroy(&fd);
141: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2],FILE_MODE_READ,&fd);
142: VecSetFromOptions(b);
143: VecLoad(b,fd);
144: }
145: } else { /* rhs is stored in the same file as matrix */
146: VecSetFromOptions(b);
147: VecLoad(b,fd);
148: }
149: PetscViewerDestroy(&fd);
151: /* Make A singular for testing zero-pivot of ilu factorization */
152: /* Example: ./ex10 -f0 <datafile> -test_zeropivot -set_row_zero -pc_factor_shift_type <shift_type> */
153: flg = PETSC_FALSE;
154: PetscOptionsGetBool(NULL,NULL, "-test_zeropivot", &flg,NULL);
155: if (flg) {
156: PetscInt row=0;
157: PetscBool flg1=PETSC_FALSE;
158: PetscScalar zero=0.0;
160: PetscOptionsGetBool(NULL,NULL, "-set_row_zero", &flg1,NULL);
161: if (flg1) { /* set a row as zeros */
162: MatSetOption(A,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE);
163: MatZeroRows(A,1,&row,0.0,NULL,NULL);
164: } else if (!rank) {
165: /* set (row,row) entry as zero */
166: MatSetValues(A,1,&row,1,&row,&zero,INSERT_VALUES);
167: }
168: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
169: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
170: }
172: /* Check whether A is symmetric */
173: flg = PETSC_FALSE;
174: PetscOptionsGetBool(NULL,NULL, "-check_symmetry", &flg,NULL);
175: if (flg) {
176: Mat Atrans;
177: MatTranspose(A, MAT_INITIAL_MATRIX,&Atrans);
178: MatEqual(A, Atrans, &isSymmetric);
179: if (isSymmetric) {
180: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
181: } else {
182: PetscPrintf(PETSC_COMM_WORLD,"Warning: A is non-symmetric \n");
183: }
184: MatDestroy(&Atrans);
185: }
187: /*
188: If the loaded matrix is larger than the vector (due to being padded
189: to match the block size of the system), then create a new padded vector.
190: */
192: MatGetLocalSize(A,&m,&n);
193: /* if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);*/
194: MatGetSize(A,&M,NULL);
195: VecGetSize(b,&m);
196: if (M != m) { /* Create a new vector b by padding the old one */
197: PetscInt j,mvec,start,end,indx;
198: Vec tmp;
199: PetscScalar *bold;
201: VecCreate(PETSC_COMM_WORLD,&tmp);
202: VecSetSizes(tmp,n,PETSC_DECIDE);
203: VecSetFromOptions(tmp);
204: VecGetOwnershipRange(b,&start,&end);
205: VecGetLocalSize(b,&mvec);
206: VecGetArray(b,&bold);
207: for (j=0; j<mvec; j++) {
208: indx = start+j;
209: VecSetValues(tmp,1,&indx,bold+j,INSERT_VALUES);
210: }
211: VecRestoreArray(b,&bold);
212: VecDestroy(&b);
213: VecAssemblyBegin(tmp);
214: VecAssemblyEnd(tmp);
215: b = tmp;
216: }
218: MatCreateVecs(A,&x,NULL);
219: VecDuplicate(b,&u);
220: if (initialguessfile) {
221: PetscViewer viewer2;
222: PetscViewerBinaryOpen(PETSC_COMM_WORLD,initialguessfilename,FILE_MODE_READ,&viewer2);
223: VecLoad(x,viewer2);
224: PetscViewerDestroy(&viewer2);
225: initialguess = PETSC_TRUE;
226: } else if (initialguess) {
227: VecSet(x,1.0);
228: } else {
229: VecSet(x,0.0);
230: }
233: /* Check scaling in A */
234: flg = PETSC_FALSE;
235: PetscOptionsGetBool(NULL,NULL, "-check_scaling", &flg,NULL);
236: if (flg) {
237: Vec max, min;
238: PetscInt idx;
239: PetscReal val;
241: VecDuplicate(x, &max);
242: VecDuplicate(x, &min);
243: MatGetRowMaxAbs(A, max, NULL);
244: MatGetRowMinAbs(A, min, NULL);
245: {
246: PetscViewer viewer;
248: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "max.data", &viewer);
249: VecView(max, viewer);
250: PetscViewerDestroy(&viewer);
251: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "min.data", &viewer);
252: VecView(min, viewer);
253: PetscViewerDestroy(&viewer);
254: }
255: VecView(max, PETSC_VIEWER_DRAW_WORLD);
256: VecMax(max, &idx, &val);
257: PetscPrintf(PETSC_COMM_WORLD, "Largest max row element %g at row %D\n", (double)val, idx);
258: VecView(min, PETSC_VIEWER_DRAW_WORLD);
259: VecMin(min, &idx, &val);
260: PetscPrintf(PETSC_COMM_WORLD, "Smallest min row element %g at row %D\n", (double)val, idx);
261: VecMin(max, &idx, &val);
262: PetscPrintf(PETSC_COMM_WORLD, "Smallest max row element %g at row %D\n", (double)val, idx);
263: VecPointwiseDivide(max, max, min);
264: VecMax(max, &idx, &val);
265: PetscPrintf(PETSC_COMM_WORLD, "Largest row ratio %g at row %D\n", (double)val, idx);
266: VecView(max, PETSC_VIEWER_DRAW_WORLD);
267: VecDestroy(&max);
268: VecDestroy(&min);
269: }
271: /* MatView(A,PETSC_VIEWER_STDOUT_WORLD); */
272: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
273: Setup solve for system
274: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
275: /*
276: Conclude profiling last stage; begin profiling next stage.
277: */
278: PetscPreLoadStage("KSPSetUpSolve");
280: /*
281: Create linear solver; set operators; set runtime options.
282: */
283: KSPCreate(PETSC_COMM_WORLD,&ksp);
284: KSPSetInitialGuessNonzero(ksp,initialguess);
285: num_numfac = 1;
286: PetscOptionsGetInt(NULL,NULL,"-num_numfac",&num_numfac,NULL);
287: while (num_numfac--) {
288: PetscBool lsqr;
289: char str[32];
290: PetscOptionsGetString(NULL,NULL,"-ksp_type",str,32,&lsqr);
291: if (lsqr) {
292: PetscStrcmp("lsqr",str,&lsqr);
293: }
294: if (lsqr) {
295: Mat BtB;
296: MatTransposeMatMult(A,A,MAT_INITIAL_MATRIX,4,&BtB);
297: KSPSetOperators(ksp,A,BtB);
298: MatDestroy(&BtB);
299: } else {
300: KSPSetOperators(ksp,A,A);
301: }
302: KSPSetFromOptions(ksp);
304: /*
305: Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
306: enable more precise profiling of setting up the preconditioner.
307: These calls are optional, since both will be called within
308: KSPSolve() if they haven't been called already.
309: */
310: KSPSetUp(ksp);
311: KSPSetUpOnBlocks(ksp);
313: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
314: Solve system
315: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
317: /*
318: Solve linear system;
319: */
320: if (trans) {
321: KSPSolveTranspose(ksp,b,x);
322: KSPGetIterationNumber(ksp,&its);
323: } else {
324: PetscInt num_rhs=1;
325: PetscOptionsGetInt(NULL,NULL,"-num_rhs",&num_rhs,NULL);
326: cknorm = PETSC_FALSE;
327: PetscOptionsGetBool(NULL,NULL,"-cknorm",&cknorm,NULL);
328: while (num_rhs--) {
329: if (num_rhs == 1) VecSet(x,0.0);
330: KSPSolve(ksp,b,x);
331: }
332: KSPGetIterationNumber(ksp,&its);
333: if (cknorm) { /* Check error for each rhs */
334: if (trans) {
335: MatMultTranspose(A,x,u);
336: } else {
337: MatMult(A,x,u);
338: }
339: VecAXPY(u,-1.0,b);
340: VecNorm(u,NORM_2,&norm);
341: PetscPrintf(PETSC_COMM_WORLD," Number of iterations = %3D\n",its);
342: if (!PetscIsNanScalar(norm)) {
343: if (norm < 1.e-12) {
344: PetscPrintf(PETSC_COMM_WORLD," Residual norm < 1.e-12\n");
345: } else {
346: PetscPrintf(PETSC_COMM_WORLD," Residual norm %g\n",(double)norm);
347: }
348: }
349: }
350: } /* while (num_rhs--) */
352: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
353: Check error, print output, free data structures.
354: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
356: /*
357: Check error
358: */
359: if (trans) {
360: MatMultTranspose(A,x,u);
361: } else {
362: MatMult(A,x,u);
363: }
364: VecAXPY(u,-1.0,b);
365: VecNorm(u,NORM_2,&norm);
366: /*
367: Write output (optinally using table for solver details).
368: - PetscPrintf() handles output for multiprocessor jobs
369: by printing from only one processor in the communicator.
370: - KSPView() prints information about the linear solver.
371: */
372: if (table) {
373: char *matrixname,kspinfo[120];
374: PetscViewer viewer;
376: /*
377: Open a string viewer; then write info to it.
378: */
379: PetscViewerStringOpen(PETSC_COMM_WORLD,kspinfo,120,&viewer);
380: KSPView(ksp,viewer);
381: PetscStrrchr(file[PetscPreLoadIt],'/',&matrixname);
382: PetscPrintf(PETSC_COMM_WORLD,"%-8.8s %3D %2.0e %s \n",matrixname,its,norm,kspinfo);
384: /*
385: Destroy the viewer
386: */
387: PetscViewerDestroy(&viewer);
388: } else {
389: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
390: if (!PetscIsNanScalar(norm)) {
391: if (norm < 1.e-12 && !PetscIsNanScalar((PetscScalar)norm)) {
392: PetscPrintf(PETSC_COMM_WORLD," Residual norm < 1.e-12\n");
393: } else {
394: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);
395: }
396: }
397: }
398: PetscOptionsGetString(NULL,NULL,"-solution",file[3],PETSC_MAX_PATH_LEN,&flg);
399: if (flg) {
400: PetscViewer viewer;
401: Vec xstar;
402: PetscReal norm;
404: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[3],FILE_MODE_READ,&viewer);
405: VecCreate(PETSC_COMM_WORLD,&xstar);
406: VecLoad(xstar,viewer);
407: VecAXPY(xstar, -1.0, x);
408: VecNorm(xstar, NORM_2, &norm);
409: PetscPrintf(PETSC_COMM_WORLD, "Error norm %g\n", (double)norm);
410: VecDestroy(&xstar);
411: PetscViewerDestroy(&viewer);
412: }
413: if (outputSoln) {
414: PetscViewer viewer;
416: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"solution.petsc",FILE_MODE_WRITE,&viewer);
417: VecView(x, viewer);
418: PetscViewerDestroy(&viewer);
419: }
421: flg = PETSC_FALSE;
422: PetscOptionsGetBool(NULL,NULL, "-ksp_reason", &flg,NULL);
423: if (flg) {
424: KSPConvergedReason reason;
425: KSPGetConvergedReason(ksp,&reason);
426: PetscPrintf(PETSC_COMM_WORLD,"KSPConvergedReason: %D\n", reason);
427: }
429: } /* while (num_numfac--) */
431: /*
432: Free work space. All PETSc objects should be destroyed when they
433: are no longer needed.
434: */
435: MatDestroy(&A); VecDestroy(&b);
436: VecDestroy(&u); VecDestroy(&x);
437: KSPDestroy(&ksp);
438: PetscPreLoadEnd();
439: /* -----------------------------------------------------------
440: End of linear solver loop
441: ----------------------------------------------------------- */
443: PetscFinalize();
444: return 0;
445: }