Actual source code: ex10.c
petsc-3.5.4 2015-05-23
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,"-table",&table,NULL);
60: PetscOptionsGetBool(NULL,"-trans",&trans,NULL);
61: PetscOptionsGetBool(NULL,"-initialguess",&initialguess,NULL);
62: PetscOptionsGetBool(NULL,"-output_solution",&outputSoln,NULL);
63: PetscOptionsGetString(NULL,"-initialguessfilename",initialguessfilename,PETSC_MAX_PATH_LEN,&initialguessfile);
64: PetscOptionsGetInt(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,"-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,"-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,"-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: PetscMalloc(nearnulldim*sizeof(nullvecs[0]),&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,"-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_nonzero */
153: flg = PETSC_FALSE;
154: PetscOptionsGetBool(NULL, "-test_zeropivot", &flg,NULL);
155: if (flg) {
156: PetscInt row,ncols;
157: const PetscInt *cols;
158: const PetscScalar *vals;
159: PetscBool flg1=PETSC_FALSE;
160: PetscScalar *zeros;
161: row = 0;
162: MatGetRow(A,row,&ncols,&cols,&vals);
163: PetscMalloc(sizeof(PetscScalar)*(ncols+1),&zeros);
164: PetscMemzero(zeros,(ncols+1)*sizeof(PetscScalar));
165: PetscOptionsGetBool(NULL, "-set_row_zero", &flg1,NULL);
166: if (flg1) { /* set entire row as zero */
167: MatSetValues(A,1,&row,ncols,cols,zeros,INSERT_VALUES);
168: } else { /* only set (row,row) entry as zero */
169: MatSetValues(A,1,&row,1,&row,zeros,INSERT_VALUES);
170: }
171: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
172: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
173: }
175: /* Check whether A is symmetric */
176: flg = PETSC_FALSE;
177: PetscOptionsGetBool(NULL, "-check_symmetry", &flg,NULL);
178: if (flg) {
179: Mat Atrans;
180: MatTranspose(A, MAT_INITIAL_MATRIX,&Atrans);
181: MatEqual(A, Atrans, &isSymmetric);
182: if (isSymmetric) {
183: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
184: } else {
185: PetscPrintf(PETSC_COMM_WORLD,"Warning: A is non-symmetric \n");
186: }
187: MatDestroy(&Atrans);
188: }
190: /*
191: If the loaded matrix is larger than the vector (due to being padded
192: to match the block size of the system), then create a new padded vector.
193: */
195: MatGetLocalSize(A,&m,&n);
196: /* if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);*/
197: MatGetSize(A,&M,NULL);
198: VecGetSize(b,&m);
199: if (M != m) { /* Create a new vector b by padding the old one */
200: PetscInt j,mvec,start,end,indx;
201: Vec tmp;
202: PetscScalar *bold;
204: VecCreate(PETSC_COMM_WORLD,&tmp);
205: VecSetSizes(tmp,n,PETSC_DECIDE);
206: VecSetFromOptions(tmp);
207: VecGetOwnershipRange(b,&start,&end);
208: VecGetLocalSize(b,&mvec);
209: VecGetArray(b,&bold);
210: for (j=0; j<mvec; j++) {
211: indx = start+j;
212: VecSetValues(tmp,1,&indx,bold+j,INSERT_VALUES);
213: }
214: VecRestoreArray(b,&bold);
215: VecDestroy(&b);
216: VecAssemblyBegin(tmp);
217: VecAssemblyEnd(tmp);
218: b = tmp;
219: }
221: MatGetVecs(A,&x,NULL);
222: VecDuplicate(b,&u);
223: if (initialguessfile) {
224: PetscViewer viewer2;
225: PetscViewerBinaryOpen(PETSC_COMM_WORLD,initialguessfilename,FILE_MODE_READ,&viewer2);
226: VecLoad(x,viewer2);
227: PetscViewerDestroy(&viewer2);
228: initialguess = PETSC_TRUE;
229: } else if (initialguess) {
230: VecSet(x,1.0);
231: } else {
232: VecSet(x,0.0);
233: }
236: /* Check scaling in A */
237: flg = PETSC_FALSE;
238: PetscOptionsGetBool(NULL, "-check_scaling", &flg,NULL);
239: if (flg) {
240: Vec max, min;
241: PetscInt idx;
242: PetscReal val;
244: VecDuplicate(x, &max);
245: VecDuplicate(x, &min);
246: MatGetRowMaxAbs(A, max, NULL);
247: MatGetRowMinAbs(A, min, NULL);
248: {
249: PetscViewer viewer;
251: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "max.data", &viewer);
252: VecView(max, viewer);
253: PetscViewerDestroy(&viewer);
254: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "min.data", &viewer);
255: VecView(min, viewer);
256: PetscViewerDestroy(&viewer);
257: }
258: VecView(max, PETSC_VIEWER_DRAW_WORLD);
259: VecMax(max, &idx, &val);
260: PetscPrintf(PETSC_COMM_WORLD, "Largest max row element %g at row %D\n", (double)val, idx);
261: VecView(min, PETSC_VIEWER_DRAW_WORLD);
262: VecMin(min, &idx, &val);
263: PetscPrintf(PETSC_COMM_WORLD, "Smallest min row element %g at row %D\n", (double)val, idx);
264: VecMin(max, &idx, &val);
265: PetscPrintf(PETSC_COMM_WORLD, "Smallest max row element %g at row %D\n", (double)val, idx);
266: VecPointwiseDivide(max, max, min);
267: VecMax(max, &idx, &val);
268: PetscPrintf(PETSC_COMM_WORLD, "Largest row ratio %g at row %D\n", (double)val, idx);
269: VecView(max, PETSC_VIEWER_DRAW_WORLD);
270: VecDestroy(&max);
271: VecDestroy(&min);
272: }
274: /* MatView(A,PETSC_VIEWER_STDOUT_WORLD); */
275: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
276: Setup solve for system
277: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
278: /*
279: Conclude profiling last stage; begin profiling next stage.
280: */
281: PetscPreLoadStage("KSPSetUpSolve");
283: /*
284: Create linear solver; set operators; set runtime options.
285: */
286: KSPCreate(PETSC_COMM_WORLD,&ksp);
287: KSPSetInitialGuessNonzero(ksp,initialguess);
288: num_numfac = 1;
289: PetscOptionsGetInt(NULL,"-num_numfac",&num_numfac,NULL);
290: while (num_numfac--) {
291: PetscBool lsqr;
292: char str[32];
293: PetscOptionsGetString(NULL,"-ksp_type",str,32,&lsqr);
294: if (lsqr) {
295: PetscStrcmp("lsqr",str,&lsqr);
296: }
297: if (lsqr) {
298: Mat BtB;
299: MatTransposeMatMult(A,A,MAT_INITIAL_MATRIX,4,&BtB);
300: KSPSetOperators(ksp,A,BtB);
301: MatDestroy(&BtB);
302: } else {
303: KSPSetOperators(ksp,A,A);
304: }
305: KSPSetFromOptions(ksp);
307: /*
308: Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
309: enable more precise profiling of setting up the preconditioner.
310: These calls are optional, since both will be called within
311: KSPSolve() if they haven't been called already.
312: */
313: KSPSetUp(ksp);
314: KSPSetUpOnBlocks(ksp);
316: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
317: Solve system
318: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
320: /*
321: Solve linear system;
322: */
323: if (trans) {
324: KSPSolveTranspose(ksp,b,x);
325: KSPGetIterationNumber(ksp,&its);
326: } else {
327: PetscInt num_rhs=1;
328: PetscOptionsGetInt(NULL,"-num_rhs",&num_rhs,NULL);
329: cknorm = PETSC_FALSE;
330: PetscOptionsGetBool(NULL,"-cknorm",&cknorm,NULL);
331: while (num_rhs--) {
332: if (num_rhs == 1) VecSet(x,0.0);
333: KSPSolve(ksp,b,x);
334: }
335: KSPGetIterationNumber(ksp,&its);
336: if (cknorm) { /* Check error for each rhs */
337: if (trans) {
338: MatMultTranspose(A,x,u);
339: } else {
340: MatMult(A,x,u);
341: }
342: VecAXPY(u,-1.0,b);
343: VecNorm(u,NORM_2,&norm);
344: PetscPrintf(PETSC_COMM_WORLD," Number of iterations = %3D\n",its);
345: if (norm < 1.e-12) {
346: PetscPrintf(PETSC_COMM_WORLD," Residual norm < 1.e-12\n");
347: } else {
348: PetscPrintf(PETSC_COMM_WORLD," Residual norm %g\n",(double)norm);
349: }
350: }
351: } /* while (num_rhs--) */
353: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
354: Check error, print output, free data structures.
355: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
357: /*
358: Check error
359: */
360: if (trans) {
361: MatMultTranspose(A,x,u);
362: } else {
363: MatMult(A,x,u);
364: }
365: VecAXPY(u,-1.0,b);
366: VecNorm(u,NORM_2,&norm);
367: /*
368: Write output (optinally using table for solver details).
369: - PetscPrintf() handles output for multiprocessor jobs
370: by printing from only one processor in the communicator.
371: - KSPView() prints information about the linear solver.
372: */
373: if (table) {
374: char *matrixname,kspinfo[120];
375: PetscViewer viewer;
377: /*
378: Open a string viewer; then write info to it.
379: */
380: PetscViewerStringOpen(PETSC_COMM_WORLD,kspinfo,120,&viewer);
381: KSPView(ksp,viewer);
382: PetscStrrchr(file[PetscPreLoadIt],'/',&matrixname);
383: PetscPrintf(PETSC_COMM_WORLD,"%-8.8s %3D %2.0e %s \n",matrixname,its,norm,kspinfo);
385: /*
386: Destroy the viewer
387: */
388: PetscViewerDestroy(&viewer);
389: } else {
390: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
391: if (norm < 1.e-12) {
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: PetscOptionsGetString(NULL,"-solution",file[3],PETSC_MAX_PATH_LEN,&flg);
398: if (flg) {
399: PetscViewer viewer;
400: Vec xstar;
401: PetscReal norm;
403: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[3],FILE_MODE_READ,&viewer);
404: VecCreate(PETSC_COMM_WORLD,&xstar);
405: VecLoad(xstar,viewer);
406: VecAXPY(xstar, -1.0, x);
407: VecNorm(xstar, NORM_2, &norm);
408: PetscPrintf(PETSC_COMM_WORLD, "Error norm %g\n", (double)norm);
409: VecDestroy(&xstar);
410: PetscViewerDestroy(&viewer);
411: }
412: if (outputSoln) {
413: PetscViewer viewer;
415: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"solution.petsc",FILE_MODE_WRITE,&viewer);
416: VecView(x, viewer);
417: PetscViewerDestroy(&viewer);
418: }
420: flg = PETSC_FALSE;
421: PetscOptionsGetBool(NULL, "-ksp_reason", &flg,NULL);
422: if (flg) {
423: KSPConvergedReason reason;
424: KSPGetConvergedReason(ksp,&reason);
425: PetscPrintf(PETSC_COMM_WORLD,"KSPConvergedReason: %D\n", reason);
426: }
428: } /* while (num_numfac--) */
430: /*
431: Free work space. All PETSc objects should be destroyed when they
432: are no longer needed.
433: */
434: MatDestroy(&A); VecDestroy(&b);
435: VecDestroy(&u); VecDestroy(&x);
436: KSPDestroy(&ksp);
437: PetscPreLoadEnd();
438: /* -----------------------------------------------------------
439: End of linear solver loop
440: ----------------------------------------------------------- */
442: PetscFinalize();
443: return 0;
444: }