Actual source code: ex10.c
petsc-3.4.5 2014-06-29
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: -trans : solve transpose system instead\n\n";
12: /*
13: This code can be used to test PETSc interface to other packages.\n\
14: Examples of command line options: \n\
15: ./ex10 -f0 <datafile> -ksp_type preonly \n\
16: -help -ksp_view \n\
17: -num_numfac <num_numfac> -num_rhs <num_rhs> \n\
18: -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package superlu or superlu_dist or mumps \n\
19: -ksp_type preonly -pc_type cholesky -pc_factor_mat_solver_package mumps \n\
20: mpiexec -n <np> ./ex10 -f0 <datafile> -ksp_type cg -pc_type asm -pc_asm_type basic -sub_pc_type icc -mat_type sbaij
21: \n\n";
22: */
23: /*T
24: Concepts: KSP^solving a linear system
25: Processors: n
26: T*/
28: /*
29: Include "petscksp.h" so that we can use KSP solvers. Note that this file
30: automatically includes:
31: petscsys.h - base PETSc routines petscvec.h - vectors
32: petscmat.h - matrices
33: petscis.h - index sets petscksp.h - Krylov subspace methods
34: petscviewer.h - viewers petscpc.h - preconditioners
35: */
36: #include <petscksp.h>
37: #include <petsctime.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;
52: PetscReal norm;
53: PetscLogDouble tsetup,tsetup1,tsetup2,tsolve,tsolve1,tsolve2;
54: PetscBool preload=PETSC_TRUE,isSymmetric,cknorm=PETSC_FALSE,initialguessfile = PETSC_FALSE;
55: PetscMPIInt rank;
56: char initialguessfilename[PETSC_MAX_PATH_LEN];
58: PetscInitialize(&argc,&args,(char*)0,help);
59: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
60: PetscOptionsGetBool(NULL,"-table",&table,NULL);
61: PetscOptionsGetBool(NULL,"-trans",&trans,NULL);
62: PetscOptionsGetBool(NULL,"-initialguess",&initialguess,NULL);
63: PetscOptionsGetBool(NULL,"-output_solution",&outputSoln,NULL);
64: PetscOptionsGetString(NULL,"-initialguessfilename",initialguessfilename,PETSC_MAX_PATH_LEN,&initialguessfile);
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);
112: flg = PETSC_FALSE;
113: PetscOptionsGetString(NULL,"-rhs",file[2],PETSC_MAX_PATH_LEN,&flg);
114: VecCreate(PETSC_COMM_WORLD,&b);
115: if (flg) { /* rhs is stored in a separate file */
116: if (file[2][0] == '0') {
117: PetscInt m;
118: PetscScalar one = 1.0;
119: PetscInfo(0,"Using vector of ones for RHS\n");
120: MatGetLocalSize(A,&m,NULL);
121: VecSetSizes(b,m,PETSC_DECIDE);
122: VecSetFromOptions(b);
123: VecSet(b,one);
124: } else {
125: PetscViewerDestroy(&fd);
126: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2],FILE_MODE_READ,&fd);
127: VecSetFromOptions(b);
128: VecLoad(b,fd);
129: }
130: } else { /* rhs is stored in the same file as matrix */
131: VecSetFromOptions(b);
132: VecLoad(b,fd);
133: }
134: PetscViewerDestroy(&fd);
136: /* Make A singular for testing zero-pivot of ilu factorization */
137: /* Example: ./ex10 -f0 <datafile> -test_zeropivot -set_row_zero -pc_factor_shift_nonzero */
138: flg = PETSC_FALSE;
139: PetscOptionsGetBool(NULL, "-test_zeropivot", &flg,NULL);
140: if (flg) {
141: PetscInt row,ncols;
142: const PetscInt *cols;
143: const PetscScalar *vals;
144: PetscBool flg1=PETSC_FALSE;
145: PetscScalar *zeros;
146: row = 0;
147: MatGetRow(A,row,&ncols,&cols,&vals);
148: PetscMalloc(sizeof(PetscScalar)*(ncols+1),&zeros);
149: PetscMemzero(zeros,(ncols+1)*sizeof(PetscScalar));
150: PetscOptionsGetBool(NULL, "-set_row_zero", &flg1,NULL);
151: if (flg1) { /* set entire row as zero */
152: MatSetValues(A,1,&row,ncols,cols,zeros,INSERT_VALUES);
153: } else { /* only set (row,row) entry as zero */
154: MatSetValues(A,1,&row,1,&row,zeros,INSERT_VALUES);
155: }
156: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
157: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
158: }
160: /* Check whether A is symmetric */
161: flg = PETSC_FALSE;
162: PetscOptionsGetBool(NULL, "-check_symmetry", &flg,NULL);
163: if (flg) {
164: Mat Atrans;
165: MatTranspose(A, MAT_INITIAL_MATRIX,&Atrans);
166: MatEqual(A, Atrans, &isSymmetric);
167: if (isSymmetric) {
168: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
169: } else {
170: PetscPrintf(PETSC_COMM_WORLD,"Warning: A is non-symmetric \n");
171: }
172: MatDestroy(&Atrans);
173: }
175: /*
176: If the loaded matrix is larger than the vector (due to being padded
177: to match the block size of the system), then create a new padded vector.
178: */
180: MatGetLocalSize(A,&m,&n);
181: /* if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);*/
182: MatGetSize(A,&M,NULL);
183: VecGetSize(b,&m);
184: if (M != m) { /* Create a new vector b by padding the old one */
185: PetscInt j,mvec,start,end,indx;
186: Vec tmp;
187: PetscScalar *bold;
189: VecCreate(PETSC_COMM_WORLD,&tmp);
190: VecSetSizes(tmp,n,PETSC_DECIDE);
191: VecSetFromOptions(tmp);
192: VecGetOwnershipRange(b,&start,&end);
193: VecGetLocalSize(b,&mvec);
194: VecGetArray(b,&bold);
195: for (j=0; j<mvec; j++) {
196: indx = start+j;
197: VecSetValues(tmp,1,&indx,bold+j,INSERT_VALUES);
198: }
199: VecRestoreArray(b,&bold);
200: VecDestroy(&b);
201: VecAssemblyBegin(tmp);
202: VecAssemblyEnd(tmp);
203: b = tmp;
204: }
206: MatGetVecs(A,&x,NULL);
207: VecDuplicate(b,&u);
208: if (initialguessfile) {
209: PetscViewer viewer2;
210: PetscViewerBinaryOpen(PETSC_COMM_WORLD,initialguessfilename,FILE_MODE_READ,&viewer2);
211: VecLoad(x,viewer2);
212: PetscViewerDestroy(&viewer2);
213: initialguess = PETSC_TRUE;
214: } else if (initialguess) {
215: VecSet(x,1.0);
216: } else {
217: VecSet(x,0.0);
218: }
221: /* Check scaling in A */
222: flg = PETSC_FALSE;
223: PetscOptionsGetBool(NULL, "-check_scaling", &flg,NULL);
224: if (flg) {
225: Vec max, min;
226: PetscInt idx;
227: PetscReal val;
229: VecDuplicate(x, &max);
230: VecDuplicate(x, &min);
231: MatGetRowMaxAbs(A, max, NULL);
232: MatGetRowMinAbs(A, min, NULL);
233: {
234: PetscViewer viewer;
236: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "max.data", &viewer);
237: VecView(max, viewer);
238: PetscViewerDestroy(&viewer);
239: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "min.data", &viewer);
240: VecView(min, viewer);
241: PetscViewerDestroy(&viewer);
242: }
243: VecView(max, PETSC_VIEWER_DRAW_WORLD);
244: VecMax(max, &idx, &val);
245: PetscPrintf(PETSC_COMM_WORLD, "Largest max row element %G at row %d\n", val, idx);
246: VecView(min, PETSC_VIEWER_DRAW_WORLD);
247: VecMin(min, &idx, &val);
248: PetscPrintf(PETSC_COMM_WORLD, "Smallest min row element %G at row %d\n", val, idx);
249: VecMin(max, &idx, &val);
250: PetscPrintf(PETSC_COMM_WORLD, "Smallest max row element %G at row %d\n", val, idx);
251: VecPointwiseDivide(max, max, min);
252: VecMax(max, &idx, &val);
253: PetscPrintf(PETSC_COMM_WORLD, "Largest row ratio %G at row %d\n", val, idx);
254: VecView(max, PETSC_VIEWER_DRAW_WORLD);
255: VecDestroy(&max);
256: VecDestroy(&min);
257: }
259: /* MatView(A,PETSC_VIEWER_STDOUT_WORLD); */
260: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
261: Setup solve for system
262: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
263: /*
264: Conclude profiling last stage; begin profiling next stage.
265: */
266: PetscPreLoadStage("KSPSetUpSolve");
268: /*
269: We also explicitly time this stage via PetscTime()
270: */
271: PetscTime(&tsetup1);
273: /*
274: Create linear solver; set operators; set runtime options.
275: */
276: KSPCreate(PETSC_COMM_WORLD,&ksp);
277: KSPSetInitialGuessNonzero(ksp,initialguess);
278: num_numfac = 1;
279: PetscOptionsGetInt(NULL,"-num_numfac",&num_numfac,NULL);
280: while (num_numfac--) {
281: PetscBool lsqr;
282: char str[32];
283: PetscOptionsGetString(NULL,"-ksp_type",str,32,&lsqr);
284: if (lsqr) {
285: PetscStrcmp("lsqr",str,&lsqr);
286: }
287: if (lsqr) {
288: Mat BtB;
289: MatTransposeMatMult(A,A,MAT_INITIAL_MATRIX,4,&BtB);
290: KSPSetOperators(ksp,A,BtB,SAME_NONZERO_PATTERN);
291: MatDestroy(&BtB);
292: } else {
293: KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);
294: }
295: KSPSetFromOptions(ksp);
297: /*
298: Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
299: enable more precise profiling of setting up the preconditioner.
300: These calls are optional, since both will be called within
301: KSPSolve() if they haven't been called already.
302: */
303: KSPSetUp(ksp);
304: KSPSetUpOnBlocks(ksp);
305: PetscTime(&tsetup2);
306: tsetup = tsetup2 - tsetup1;
308: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
309: Solve system
310: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
312: /*
313: Solve linear system; we also explicitly time this stage.
314: */
315: PetscTime(&tsolve1);
316: if (trans) {
317: KSPSolveTranspose(ksp,b,x);
318: KSPGetIterationNumber(ksp,&its);
319: } else {
320: PetscInt num_rhs=1;
321: PetscOptionsGetInt(NULL,"-num_rhs",&num_rhs,NULL);
322: cknorm = PETSC_FALSE;
323: PetscOptionsGetBool(NULL,"-cknorm",&cknorm,NULL);
324: while (num_rhs--) {
325: if (num_rhs == 1) VecSet(x,0.0);
326: KSPSolve(ksp,b,x);
327: }
328: KSPGetIterationNumber(ksp,&its);
329: if (cknorm) { /* Check error for each rhs */
330: if (trans) {
331: MatMultTranspose(A,x,u);
332: } else {
333: MatMult(A,x,u);
334: }
335: VecAXPY(u,-1.0,b);
336: VecNorm(u,NORM_2,&norm);
337: PetscPrintf(PETSC_COMM_WORLD," Number of iterations = %3D\n",its);
338: if (norm < 1.e-12) {
339: PetscPrintf(PETSC_COMM_WORLD," Residual norm < 1.e-12\n");
340: } else {
341: PetscPrintf(PETSC_COMM_WORLD," Residual norm %G\n",norm);
342: }
343: }
344: } /* while (num_rhs--) */
345: PetscTime(&tsolve2);
346: tsolve = tsolve2 - tsolve1;
348: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
349: Check error, print output, free data structures.
350: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
352: /*
353: Check error
354: */
355: if (trans) {
356: MatMultTranspose(A,x,u);
357: } else {
358: MatMult(A,x,u);
359: }
360: VecAXPY(u,-1.0,b);
361: VecNorm(u,NORM_2,&norm);
362: /*
363: Write output (optinally using table for solver details).
364: - PetscPrintf() handles output for multiprocessor jobs
365: by printing from only one processor in the communicator.
366: - KSPView() prints information about the linear solver.
367: */
368: if (table) {
369: char *matrixname,kspinfo[120];
370: PetscViewer viewer;
372: /*
373: Open a string viewer; then write info to it.
374: */
375: PetscViewerStringOpen(PETSC_COMM_WORLD,kspinfo,120,&viewer);
376: KSPView(ksp,viewer);
377: PetscStrrchr(file[PetscPreLoadIt],'/',&matrixname);
378: PetscPrintf(PETSC_COMM_WORLD,"%-8.8s %3D %2.0e %2.1e %2.1e %2.1e %s \n",
379: matrixname,its,norm,tsetup+tsolve,tsetup,tsolve,kspinfo);
381: /*
382: Destroy the viewer
383: */
384: PetscViewerDestroy(&viewer);
385: } else {
386: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
387: if (norm < 1.e-12) {
388: PetscPrintf(PETSC_COMM_WORLD," Residual norm < 1.e-12\n");
389: } else {
390: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %G\n",norm);
391: }
392: }
393: PetscOptionsGetString(NULL,"-solution",file[3],PETSC_MAX_PATH_LEN,&flg);
394: if (flg) {
395: PetscViewer viewer;
396: Vec xstar;
397: PetscReal norm;
399: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[3],FILE_MODE_READ,&viewer);
400: VecCreate(PETSC_COMM_WORLD,&xstar);
401: VecLoad(xstar,viewer);
402: VecAXPY(xstar, -1.0, x);
403: VecNorm(xstar, NORM_2, &norm);
404: PetscPrintf(PETSC_COMM_WORLD, "Error norm %G\n", norm);
405: VecDestroy(&xstar);
406: PetscViewerDestroy(&viewer);
407: }
408: if (outputSoln) {
409: PetscViewer viewer;
411: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"solution.petsc",FILE_MODE_WRITE,&viewer);
412: VecView(x, viewer);
413: PetscViewerDestroy(&viewer);
414: }
416: flg = PETSC_FALSE;
417: PetscOptionsGetBool(NULL, "-ksp_reason", &flg,NULL);
418: if (flg) {
419: KSPConvergedReason reason;
420: KSPGetConvergedReason(ksp,&reason);
421: PetscPrintf(PETSC_COMM_WORLD,"KSPConvergedReason: %D\n", reason);
422: }
424: } /* while (num_numfac--) */
426: /*
427: Free work space. All PETSc objects should be destroyed when they
428: are no longer needed.
429: */
430: MatDestroy(&A); VecDestroy(&b);
431: VecDestroy(&u); VecDestroy(&x);
432: KSPDestroy(&ksp);
433: PetscPreLoadEnd();
434: /* -----------------------------------------------------------
435: End of linear solver loop
436: ----------------------------------------------------------- */
438: PetscFinalize();
439: return 0;
440: }