Actual source code: ex10.c
petsc-3.3-p7 2013-05-11
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 spooles or superlu or superlu_dist or mumps \n\
19: -ksp_type preonly -pc_type cholesky -pc_factor_mat_solver_package spooles or 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>
40: int main(int argc,char **args)
41: {
42: KSP ksp; /* linear solver context */
43: Mat A; /* matrix */
44: Vec x,b,u; /* approx solution, RHS, exact solution */
45: PetscViewer fd; /* viewer */
46: char file[4][PETSC_MAX_PATH_LEN]; /* input file name */
47: PetscBool table=PETSC_FALSE,flg,trans=PETSC_FALSE,initialguess = PETSC_FALSE;
48: PetscBool outputSoln=PETSC_FALSE;
50: PetscInt its,num_numfac,m,n,M;
51: PetscReal norm;
52: PetscLogDouble tsetup,tsetup1,tsetup2,tsolve,tsolve1,tsolve2;
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(PETSC_NULL,"-table",&table,PETSC_NULL);
60: PetscOptionsGetBool(PETSC_NULL,"-trans",&trans,PETSC_NULL);
61: PetscOptionsGetBool(PETSC_NULL,"-initialguess",&initialguess,PETSC_NULL);
62: PetscOptionsGetBool(PETSC_NULL,"-output_solution",&outputSoln,PETSC_NULL);
63: PetscOptionsGetString(PETSC_NULL,"-initialguessfilename",initialguessfilename,PETSC_MAX_PATH_LEN,&initialguessfile);
65: /*
66: Determine files from which we read the two linear systems
67: (matrix and right-hand-side vector).
68: */
69: PetscOptionsGetString(PETSC_NULL,"-f",file[0],PETSC_MAX_PATH_LEN,&flg);
70: if (flg) {
71: PetscStrcpy(file[1],file[0]);
72: preload = PETSC_FALSE;
73: } else {
74: PetscOptionsGetString(PETSC_NULL,"-f0",file[0],PETSC_MAX_PATH_LEN,&flg);
75: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f0 or -f option");
76: PetscOptionsGetString(PETSC_NULL,"-f1",file[1],PETSC_MAX_PATH_LEN,&flg);
77: if (!flg) {preload = PETSC_FALSE;} /* don't bother with second system */
78: }
80: /* -----------------------------------------------------------
81: Beginning of linear solver loop
82: ----------------------------------------------------------- */
83: /*
84: Loop through the linear solve 2 times.
85: - The intention here is to preload and solve a small system;
86: then load another (larger) system and solve it as well.
87: This process preloads the instructions with the smaller
88: system so that more accurate performance monitoring (via
89: -log_summary) can be done with the larger one (that actually
90: is the system of interest).
91: */
92: PetscPreLoadBegin(preload,"Load system");
94: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
95: Load system
96: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
98: /*
99: Open binary file. Note that we use FILE_MODE_READ to indicate
100: reading from this file.
101: */
102: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[PetscPreLoadIt],FILE_MODE_READ,&fd);
103:
104: /*
105: Load the matrix and vector; then destroy the viewer.
106: */
107: MatCreate(PETSC_COMM_WORLD,&A);
108: MatSetFromOptions(A);
109: MatLoad(A,fd);
111: if (!preload){
112: flg = PETSC_FALSE;
113: PetscOptionsGetString(PETSC_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,PETSC_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 {
131: VecSetFromOptions(b);
132: VecLoad(b,fd);
133: }
134: }
135: PetscViewerDestroy(&fd);
137: /* Make A singular for testing zero-pivot of ilu factorization */
138: /* Example: ./ex10 -f0 <datafile> -test_zeropivot -set_row_zero -pc_factor_shift_nonzero */
139: flg = PETSC_FALSE;
140: PetscOptionsGetBool(PETSC_NULL, "-test_zeropivot", &flg,PETSC_NULL);
141: if (flg) {
142: PetscInt row,ncols;
143: const PetscInt *cols;
144: const PetscScalar *vals;
145: PetscBool flg1=PETSC_FALSE;
146: PetscScalar *zeros;
147: row = 0;
148: MatGetRow(A,row,&ncols,&cols,&vals);
149: PetscMalloc(sizeof(PetscScalar)*(ncols+1),&zeros);
150: PetscMemzero(zeros,(ncols+1)*sizeof(PetscScalar));
151: PetscOptionsGetBool(PETSC_NULL, "-set_row_zero", &flg1,PETSC_NULL);
152: if (flg1){ /* set entire row as zero */
153: MatSetValues(A,1,&row,ncols,cols,zeros,INSERT_VALUES);
154: } else { /* only set (row,row) entry as zero */
155: MatSetValues(A,1,&row,1,&row,zeros,INSERT_VALUES);
156: }
157: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
158: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
159: }
161: /* Check whether A is symmetric */
162: flg = PETSC_FALSE;
163: PetscOptionsGetBool(PETSC_NULL, "-check_symmetry", &flg,PETSC_NULL);
164: if (flg) {
165: Mat Atrans;
166: MatTranspose(A, MAT_INITIAL_MATRIX,&Atrans);
167: MatEqual(A, Atrans, &isSymmetric);
168: if (isSymmetric) {
169: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
170: } else {
171: PetscPrintf(PETSC_COMM_WORLD,"Warning: A is non-symmetric \n");
172: }
173: MatDestroy(&Atrans);
174: }
176: /*
177: If the loaded matrix is larger than the vector (due to being padded
178: to match the block size of the system), then create a new padded vector.
179: */
180:
181: MatGetLocalSize(A,&m,&n);
182: /* if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);*/
183: MatGetSize(A,&M,PETSC_NULL);
184: VecGetSize(b,&m);
185: if (M != m) { /* Create a new vector b by padding the old one */
186: PetscInt j,mvec,start,end,indx;
187: Vec tmp;
188: PetscScalar *bold;
190: VecCreate(PETSC_COMM_WORLD,&tmp);
191: VecSetSizes(tmp,n,PETSC_DECIDE);
192: VecSetFromOptions(tmp);
193: VecGetOwnershipRange(b,&start,&end);
194: VecGetLocalSize(b,&mvec);
195: VecGetArray(b,&bold);
196: for (j=0; j<mvec; j++) {
197: indx = start+j;
198: VecSetValues(tmp,1,&indx,bold+j,INSERT_VALUES);
199: }
200: VecRestoreArray(b,&bold);
201: VecDestroy(&b);
202: VecAssemblyBegin(tmp);
203: VecAssemblyEnd(tmp);
204: b = tmp;
205: }
206:
207: MatGetVecs(A,&x,PETSC_NULL);
208: VecDuplicate(b,&u);
209: if (initialguessfile) {
210: PetscViewer viewer2;
211: PetscViewerBinaryOpen(PETSC_COMM_WORLD,initialguessfilename,FILE_MODE_READ,&viewer2);
212: VecLoad(x,viewer2);
213: PetscViewerDestroy(&viewer2);
214: initialguess = PETSC_TRUE;
215: } else if (initialguess) {
216: VecSet(x,1.0);
217: } else {
218: VecSet(x,0.0);
219: }
222: /* Check scaling in A */
223: flg = PETSC_FALSE;
224: PetscOptionsGetBool(PETSC_NULL, "-check_scaling", &flg,PETSC_NULL);
225: if (flg) {
226: Vec max, min;
227: PetscInt idx;
228: PetscReal val;
230: VecDuplicate(x, &max);
231: VecDuplicate(x, &min);
232: MatGetRowMaxAbs(A, max, PETSC_NULL);
233: MatGetRowMinAbs(A, min, PETSC_NULL);
234: {
235: PetscViewer viewer;
237: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "max.data", &viewer);
238: VecView(max, viewer);
239: PetscViewerDestroy(&viewer);
240: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "min.data", &viewer);
241: VecView(min, viewer);
242: PetscViewerDestroy(&viewer);
243: }
244: VecView(max, PETSC_VIEWER_DRAW_WORLD);
245: VecMax(max, &idx, &val);
246: PetscPrintf(PETSC_COMM_WORLD, "Largest max row element %G at row %d\n", val, idx);
247: VecView(min, PETSC_VIEWER_DRAW_WORLD);
248: VecMin(min, &idx, &val);
249: PetscPrintf(PETSC_COMM_WORLD, "Smallest min row element %G at row %d\n", val, idx);
250: VecMin(max, &idx, &val);
251: PetscPrintf(PETSC_COMM_WORLD, "Smallest max row element %G at row %d\n", val, idx);
252: VecPointwiseDivide(max, max, min);
253: VecMax(max, &idx, &val);
254: PetscPrintf(PETSC_COMM_WORLD, "Largest row ratio %G at row %d\n", val, idx);
255: VecView(max, PETSC_VIEWER_DRAW_WORLD);
256: VecDestroy(&max);
257: VecDestroy(&min);
258: }
260: // MatView(A,PETSC_VIEWER_STDOUT_WORLD);
261: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
262: Setup solve for system
263: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
264: /*
265: Conclude profiling last stage; begin profiling next stage.
266: */
267: PetscPreLoadStage("KSPSetUpSolve");
269: /*
270: We also explicitly time this stage via PetscGetTime()
271: */
272: PetscGetTime(&tsetup1);
274: /*
275: Create linear solver; set operators; set runtime options.
276: */
277: KSPCreate(PETSC_COMM_WORLD,&ksp);
278: KSPSetInitialGuessNonzero(ksp,initialguess);
279: num_numfac = 1;
280: PetscOptionsGetInt(PETSC_NULL,"-num_numfac",&num_numfac,PETSC_NULL);
281: while ( num_numfac-- ){
282: PetscBool lsqr;
283: char str[32];
284: PetscOptionsGetString(PETSC_NULL,"-ksp_type",str,32,&lsqr);
285: if (lsqr) {
286: PetscStrcmp("lsqr",str,&lsqr);
287: }
288: if (lsqr) {
289: Mat BtB;
290: MatTransposeMatMult(A,A,MAT_INITIAL_MATRIX,4,&BtB);
291: KSPSetOperators(ksp,A,BtB,SAME_NONZERO_PATTERN);
292: MatDestroy(&BtB);
293: } else {
294: KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);
295: }
296: KSPSetFromOptions(ksp);
298: /*
299: Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
300: enable more precise profiling of setting up the preconditioner.
301: These calls are optional, since both will be called within
302: KSPSolve() if they haven't been called already.
303: */
304: KSPSetUp(ksp);
305: KSPSetUpOnBlocks(ksp);
306: PetscGetTime(&tsetup2);
307: tsetup = tsetup2 - tsetup1;
309: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
310: Solve system
311: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
313: /*
314: Solve linear system; we also explicitly time this stage.
315: */
316: PetscGetTime(&tsolve1);
317: if (trans) {
318: KSPSolveTranspose(ksp,b,x);
319: KSPGetIterationNumber(ksp,&its);
320: } else {
321: PetscInt num_rhs=1;
322: PetscOptionsGetInt(PETSC_NULL,"-num_rhs",&num_rhs,PETSC_NULL);
323: cknorm = PETSC_FALSE;
324: PetscOptionsGetBool(PETSC_NULL,"-cknorm",&cknorm,PETSC_NULL);
325: while ( num_rhs-- ) {
326: if (num_rhs == 1) VecSet(x,0.0);
327: KSPSolve(ksp,b,x);
328: }
329: KSPGetIterationNumber(ksp,&its);
330: if (cknorm){ /* Check error for each rhs */
331: if (trans) {
332: MatMultTranspose(A,x,u);
333: } else {
334: MatMult(A,x,u);
335: }
336: VecAXPY(u,-1.0,b);
337: VecNorm(u,NORM_2,&norm);
338: PetscPrintf(PETSC_COMM_WORLD," Number of iterations = %3D\n",its);
339: if (norm < 1.e-12) {
340: PetscPrintf(PETSC_COMM_WORLD," Residual norm < 1.e-12\n");
341: } else {
342: PetscPrintf(PETSC_COMM_WORLD," Residual norm %G\n",norm);
343: }
344: }
345: } /* while ( num_rhs-- ) */
346: PetscGetTime(&tsolve2);
347: tsolve = tsolve2 - tsolve1;
349: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
350: Check error, print output, free data structures.
351: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
353: /*
354: Check error
355: */
356: if (trans) {
357: MatMultTranspose(A,x,u);
358: } else {
359: MatMult(A,x,u);
360: }
361: VecAXPY(u,-1.0,b);
362: VecNorm(u,NORM_2,&norm);
363: /*
364: Write output (optinally using table for solver details).
365: - PetscPrintf() handles output for multiprocessor jobs
366: by printing from only one processor in the communicator.
367: - KSPView() prints information about the linear solver.
368: */
369: if (table) {
370: char *matrixname,kspinfo[120];
371: PetscViewer viewer;
373: /*
374: Open a string viewer; then write info to it.
375: */
376: PetscViewerStringOpen(PETSC_COMM_WORLD,kspinfo,120,&viewer);
377: KSPView(ksp,viewer);
378: PetscStrrchr(file[PetscPreLoadIt],'/',&matrixname);
379: PetscPrintf(PETSC_COMM_WORLD,"%-8.8s %3D %2.0e %2.1e %2.1e %2.1e %s \n",
380: matrixname,its,norm,tsetup+tsolve,tsetup,tsolve,kspinfo);
382: /*
383: Destroy the viewer
384: */
385: PetscViewerDestroy(&viewer);
386: } else {
387: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
388: if (norm < 1.e-12) {
389: PetscPrintf(PETSC_COMM_WORLD," Residual norm < 1.e-12\n");
390: } else {
391: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %G\n",norm);
392: }
393: }
394: PetscOptionsGetString(PETSC_NULL,"-solution",file[3],PETSC_MAX_PATH_LEN,&flg);
395: if (flg) {
396: PetscViewer viewer;
397: Vec xstar;
398: PetscReal norm;
400: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[3],FILE_MODE_READ,&viewer);
401: VecCreate(PETSC_COMM_WORLD,&xstar);
402: VecLoad(xstar,viewer);
403: VecAXPY(xstar, -1.0, x);
404: VecNorm(xstar, NORM_2, &norm);
405: PetscPrintf(PETSC_COMM_WORLD, "Error norm %G\n", norm);
406: VecDestroy(&xstar);
407: PetscViewerDestroy(&viewer);
408: }
409: if (outputSoln) {
410: PetscViewer viewer;
412: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"solution.petsc",FILE_MODE_WRITE,&viewer);
413: VecView(x, viewer);
414: PetscViewerDestroy(&viewer);
415: }
417: flg = PETSC_FALSE;
418: PetscOptionsGetBool(PETSC_NULL, "-ksp_reason", &flg,PETSC_NULL);
419: if (flg){
420: KSPConvergedReason reason;
421: KSPGetConvergedReason(ksp,&reason);
422: PetscPrintf(PETSC_COMM_WORLD,"KSPConvergedReason: %D\n", reason);
423: }
424:
425: } /* while ( num_numfac-- ) */
427: /*
428: Free work space. All PETSc objects should be destroyed when they
429: are no longer needed.
430: */
431: MatDestroy(&A); VecDestroy(&b);
432: VecDestroy(&u); VecDestroy(&x);
433: KSPDestroy(&ksp);
434: PetscPreLoadEnd();
435: /* -----------------------------------------------------------
436: End of linear solver loop
437: ----------------------------------------------------------- */
439: PetscFinalize();
440: return 0;
441: }