Actual source code: ex30.c
petsc-3.4.2 2013-07-02
2: static char help[] = "Reads a PETSc matrix and vector from a file and solves a linear system.\n\
3: It is copied and intended to move dirty codes from ksp/examples/tutorials/ex10.c and simplify ex10.c.\n\
4: Input parameters include\n\
5: -f0 <input_file> : first file to load (small system)\n\
6: -f1 <input_file> : second file to load (larger system)\n\n\
7: -trans : solve transpose system instead\n\n";
8: /*
9: This code can be used to test PETSc interface to other packages.\n\
10: Examples of command line options: \n\
11: ex30 -f0 <datafile> -ksp_type preonly \n\
12: -help -ksp_view \n\
13: -num_numfac <num_numfac> -num_rhs <num_rhs> \n\
14: -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package superlu or superlu_dist or mumps \n\
15: -ksp_type preonly -pc_type cholesky -pc_factor_mat_solver_package mumps \n\
16: mpiexec -n <np> ex30 -f0 <datafile> -ksp_type cg -pc_type asm -pc_asm_type basic -sub_pc_type icc -mat_type sbaij
18: ./ex30 -f0 $D/small -mat_sigma -3.999999999999999 -ksp_type fgmres -pc_type lu -pc_factor_mat_solver_package superlu -mat_superlu_conditionnumber -ckerror -mat_superlu_diagpivotthresh 0
19: ./ex30 -f0 $D/small -mat_sigma -3.999999999999999 -ksp_type fgmres -pc_type hypre -pc_hypre_type boomeramg -ksp_type fgmres -ckError
20: ./ex30 -f0 $D/small -mat_sigma -3.999999999999999 -ksp_type fgmres -pc_type lu -pc_factor_mat_solver_package petsc -pc_factor_shift_type NONZERO -pc_factor_shift_amount 1.e-5 -ckerror
21: \n\n";
22: */
23: /*T
24: Concepts: KSP solving a linear system
25: Processors: n
26: T*/
28: #include <petscksp.h>
29: #include <petsctime.h>
33: int main(int argc,char **args)
34: {
35: KSP ksp;
36: Mat A,B;
37: Vec x,b,u,b2; /* approx solution, RHS, exact solution */
38: PetscViewer fd; /* viewer */
39: char file[4][PETSC_MAX_PATH_LEN]; /* input file name */
40: PetscBool table = PETSC_FALSE,flg,flgB=PETSC_FALSE,trans=PETSC_FALSE,partition=PETSC_FALSE,initialguess = PETSC_FALSE;
41: PetscBool outputSoln=PETSC_FALSE;
43: PetscInt its,num_numfac,n,M;
44: PetscReal rnorm,enorm;
45: PetscLogDouble tsetup,tsetup1,tsetup2,tsolve,tsolve1,tsolve2;
46: PetscBool preload=PETSC_TRUE,diagonalscale,isSymmetric,ckrnorm=PETSC_TRUE,Test_MatDuplicate=PETSC_FALSE,ckerror=PETSC_FALSE;
47: PetscMPIInt rank;
48: PetscScalar sigma;
49: PetscInt m;
51: PetscInitialize(&argc,&args,(char*)0,help);
52: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
53: PetscOptionsGetBool(NULL,"-table",&table,NULL);
54: PetscOptionsGetBool(NULL,"-trans",&trans,NULL);
55: PetscOptionsGetBool(NULL,"-partition",&partition,NULL);
56: PetscOptionsGetBool(NULL,"-initialguess",&initialguess,NULL);
57: PetscOptionsGetBool(NULL,"-output_solution",&outputSoln,NULL);
58: PetscOptionsGetBool(NULL,"-ckrnorm",&ckrnorm,NULL);
59: PetscOptionsGetBool(NULL,"-ckerror",&ckerror,NULL);
61: /*
62: Determine files from which we read the two linear systems
63: (matrix and right-hand-side vector).
64: */
65: PetscOptionsGetString(NULL,"-f",file[0],PETSC_MAX_PATH_LEN,&flg);
66: if (flg) {
67: PetscStrcpy(file[1],file[0]);
68: preload = PETSC_FALSE;
69: } else {
70: PetscOptionsGetString(NULL,"-f0",file[0],PETSC_MAX_PATH_LEN,&flg);
71: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f0 or -f option");
72: PetscOptionsGetString(NULL,"-f1",file[1],PETSC_MAX_PATH_LEN,&flg);
73: if (!flg) preload = PETSC_FALSE; /* don't bother with second system */
74: }
76: /* -----------------------------------------------------------
77: Beginning of linear solver loop
78: ----------------------------------------------------------- */
79: /*
80: Loop through the linear solve 2 times.
81: - The intention here is to preload and solve a small system;
82: then load another (larger) system and solve it as well.
83: This process preloads the instructions with the smaller
84: system so that more accurate performance monitoring (via
85: -log_summary) can be done with the larger one (that actually
86: is the system of interest).
87: */
88: PetscPreLoadBegin(preload,"Load system");
90: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
91: Load system
92: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
94: /*
95: Open binary file. Note that we use FILE_MODE_READ to indicate
96: reading from this file.
97: */
98: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[PetscPreLoadIt],FILE_MODE_READ,&fd);
100: /*
101: Load the matrix and vector; then destroy the viewer.
102: */
103: MatCreate(PETSC_COMM_WORLD,&A);
104: MatLoad(A,fd);
106: if (!preload) {
107: flg = PETSC_FALSE;
108: PetscOptionsGetString(NULL,"-rhs",file[2],PETSC_MAX_PATH_LEN,&flg);
109: if (flg) { /* rhs is stored in a separate file */
110: PetscViewerDestroy(&fd);
111: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2],FILE_MODE_READ,&fd);
112: } else {
113: /* if file contains no RHS, then use a vector of all ones */
114: PetscInfo(0,"Using vector of ones for RHS\n");
115: MatGetLocalSize(A,&m,NULL);
116: VecCreate(PETSC_COMM_WORLD,&b);
117: VecSetSizes(b,m,PETSC_DECIDE);
118: VecSetFromOptions(b);
119: VecSet(b,1.0);
120: PetscObjectSetName((PetscObject)b, "Rhs vector");
121: }
122: }
123: PetscViewerDestroy(&fd);
125: /* Test MatDuplicate() */
126: if (Test_MatDuplicate) {
127: MatDuplicate(A,MAT_COPY_VALUES,&B);
128: MatEqual(A,B,&flg);
129: if (!flg) {
130: PetscPrintf(PETSC_COMM_WORLD," A != B \n");
131: }
132: MatDestroy(&B);
133: }
135: /* Add a shift to A */
136: PetscOptionsGetScalar(NULL,"-mat_sigma",&sigma,&flg);
137: if (flg) {
138: PetscOptionsGetString(NULL,"-fB",file[2],PETSC_MAX_PATH_LEN,&flgB);
139: if (flgB) {
140: /* load B to get A = A + sigma*B */
141: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2],FILE_MODE_READ,&fd);
142: MatCreate(PETSC_COMM_WORLD,&B);
143: MatSetOptionsPrefix(B,"B_");
144: MatLoad(B,fd);
145: PetscViewerDestroy(&fd);
146: MatAXPY(A,sigma,B,DIFFERENT_NONZERO_PATTERN); /* A <- sigma*B + A */
147: } else {
148: MatShift(A,sigma);
149: }
150: }
152: /* Make A singular for testing zero-pivot of ilu factorization */
153: /* Example: ./ex30 -f0 <datafile> -test_zeropivot -set_row_zero -pc_factor_shift_nonzero */
154: flg = PETSC_FALSE;
155: PetscOptionsGetBool(NULL, "-test_zeropivot", &flg,NULL);
156: if (flg) {
157: PetscInt row,ncols;
158: const PetscInt *cols;
159: const PetscScalar *vals;
160: PetscBool flg1=PETSC_FALSE;
161: PetscScalar *zeros;
162: row = 0;
163: MatGetRow(A,row,&ncols,&cols,&vals);
164: PetscMalloc(sizeof(PetscScalar)*(ncols+1),&zeros);
165: PetscMemzero(zeros,(ncols+1)*sizeof(PetscScalar));
166: flg1 = PETSC_FALSE;
167: PetscOptionsGetBool(NULL, "-set_row_zero", &flg1,NULL);
168: if (flg1) { /* set entire row as zero */
169: MatSetValues(A,1,&row,ncols,cols,zeros,INSERT_VALUES);
170: } else { /* only set (row,row) entry as zero */
171: MatSetValues(A,1,&row,1,&row,zeros,INSERT_VALUES);
172: }
173: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
174: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
175: }
177: /* Check whether A is symmetric */
178: flg = PETSC_FALSE;
179: PetscOptionsGetBool(NULL, "-check_symmetry", &flg,NULL);
180: if (flg) {
181: Mat Atrans;
182: MatTranspose(A, MAT_INITIAL_MATRIX,&Atrans);
183: MatEqual(A, Atrans, &isSymmetric);
184: if (isSymmetric) {
185: PetscPrintf(PETSC_COMM_WORLD,"A is symmetric \n");
186: } else {
187: PetscPrintf(PETSC_COMM_WORLD,"A is non-symmetric \n");
188: }
189: MatDestroy(&Atrans);
190: }
192: /*
193: If the loaded matrix is larger than the vector (due to being padded
194: to match the block size of the system), then create a new padded vector.
195: */
197: MatGetLocalSize(A,&m,&n);
198: if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);
199: MatGetSize(A,&M,NULL);
200: VecGetSize(b,&m);
201: if (M != m) { /* Create a new vector b by padding the old one */
202: PetscInt j,mvec,start,end,indx;
203: Vec tmp;
204: PetscScalar *bold;
206: VecCreate(PETSC_COMM_WORLD,&tmp);
207: VecSetSizes(tmp,n,PETSC_DECIDE);
208: VecSetFromOptions(tmp);
209: VecGetOwnershipRange(b,&start,&end);
210: VecGetLocalSize(b,&mvec);
211: VecGetArray(b,&bold);
212: for (j=0; j<mvec; j++) {
213: indx = start+j;
214: VecSetValues(tmp,1,&indx,bold+j,INSERT_VALUES);
215: }
216: VecRestoreArray(b,&bold);
217: VecDestroy(&b);
218: VecAssemblyBegin(tmp);
219: VecAssemblyEnd(tmp);
220: b = tmp;
221: }
222: VecDuplicate(b,&b2);
223: VecDuplicate(b,&x);
224: PetscObjectSetName((PetscObject)x, "Solution vector");
225: VecDuplicate(b,&u);
226: PetscObjectSetName((PetscObject)u, "True Solution vector");
227: VecSet(x,0.0);
229: if (ckerror) { /* Set true solution */
230: VecSet(u,1.0);
231: MatMult(A,u,b);
232: }
234: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
235: Setup solve for system
236: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
238: if (partition) {
239: MatPartitioning mpart;
240: IS mis,nis,is;
241: PetscInt *count;
242: PetscMPIInt size;
243: Mat BB;
244: MPI_Comm_size(PETSC_COMM_WORLD,&size);
245: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
246: PetscMalloc(size*sizeof(PetscInt),&count);
247: MatPartitioningCreate(PETSC_COMM_WORLD, &mpart);
248: MatPartitioningSetAdjacency(mpart, A);
249: /* MatPartitioningSetVertexWeights(mpart, weight); */
250: MatPartitioningSetFromOptions(mpart);
251: MatPartitioningApply(mpart, &mis);
252: MatPartitioningDestroy(&mpart);
253: ISPartitioningToNumbering(mis,&nis);
254: ISPartitioningCount(mis,size,count);
255: ISDestroy(&mis);
256: ISInvertPermutation(nis, count[rank], &is);
257: PetscFree(count);
258: ISDestroy(&nis);
259: ISSort(is);
260: MatGetSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&BB);
262: /* need to move the vector also */
263: ISDestroy(&is);
264: MatDestroy(&A);
265: A = BB;
266: }
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--) {
282: KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);
283: KSPSetFromOptions(ksp);
285: /*
286: Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
287: enable more precise profiling of setting up the preconditioner.
288: These calls are optional, since both will be called within
289: KSPSolve() if they haven't been called already.
290: */
291: KSPSetUp(ksp);
292: KSPSetUpOnBlocks(ksp);
293: PetscTime(&tsetup2);
294: tsetup = tsetup2 - tsetup1;
296: /*
297: Tests "diagonal-scaling of preconditioned residual norm" as used
298: by many ODE integrator codes including SUNDIALS. Note this is different
299: than diagonally scaling the matrix before computing the preconditioner
300: */
301: diagonalscale = PETSC_FALSE;
302: PetscOptionsGetBool(NULL,"-diagonal_scale",&diagonalscale,NULL);
303: if (diagonalscale) {
304: PC pc;
305: PetscInt j,start,end,n;
306: Vec scale;
308: KSPGetPC(ksp,&pc);
309: VecGetSize(x,&n);
310: VecDuplicate(x,&scale);
311: VecGetOwnershipRange(scale,&start,&end);
312: for (j=start; j<end; j++) {
313: VecSetValue(scale,j,((PetscReal)(j+1))/((PetscReal)n),INSERT_VALUES);
314: }
315: VecAssemblyBegin(scale);
316: VecAssemblyEnd(scale);
317: PCSetDiagonalScale(pc,scale);
318: VecDestroy(&scale);
319: }
321: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
322: Solve system
323: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
324: /*
325: Solve linear system; we also explicitly time this stage.
326: */
327: PetscTime(&tsolve1);
328: if (trans) {
329: KSPSolveTranspose(ksp,b,x);
330: KSPGetIterationNumber(ksp,&its);
331: } else {
332: PetscInt num_rhs=1;
333: PetscOptionsGetInt(NULL,"-num_rhs",&num_rhs,NULL);
335: while (num_rhs--) {
336: KSPSolve(ksp,b,x);
337: }
338: KSPGetIterationNumber(ksp,&its);
339: if (ckrnorm) { /* Check residual for each rhs */
340: if (trans) {
341: MatMultTranspose(A,x,b2);
342: } else {
343: MatMult(A,x,b2);
344: }
345: VecAXPY(b2,-1.0,b);
346: VecNorm(b2,NORM_2,&rnorm);
347: PetscPrintf(PETSC_COMM_WORLD," Number of iterations = %3D\n",its);
348: PetscPrintf(PETSC_COMM_WORLD," Residual norm %G\n",rnorm);
349: }
350: if (ckerror && !trans) { /* Check error for each rhs */
351: /* VecView(x,PETSC_VIEWER_STDOUT_WORLD); */
352: VecAXPY(u,-1.0,x);
353: VecNorm(u,NORM_2,&enorm);
354: PetscPrintf(PETSC_COMM_WORLD," Error norm %G\n",enorm);
355: }
357: } /* while (num_rhs--) */
358: PetscTime(&tsolve2);
359: tsolve = tsolve2 - tsolve1;
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 && ckrnorm) {
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,rnorm,tsetup+tsolve,tsetup,tsolve,kspinfo);
381: /*
382: Destroy the viewer
383: */
384: PetscViewerDestroy(&viewer);
385: }
387: PetscOptionsGetString(NULL,"-solution",file[3],PETSC_MAX_PATH_LEN,&flg);
388: if (flg) {
389: PetscViewer viewer;
390: Vec xstar;
392: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[3],FILE_MODE_READ,&viewer);
393: VecCreate(PETSC_COMM_WORLD,&xstar);
394: VecLoad(xstar,viewer);
395: VecAXPY(xstar, -1.0, x);
396: VecNorm(xstar, NORM_2, &enorm);
397: PetscPrintf(PETSC_COMM_WORLD, "Error norm %G\n", enorm);
398: VecDestroy(&xstar);
399: PetscViewerDestroy(&viewer);
400: }
401: if (outputSoln) {
402: PetscViewer viewer;
404: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"solution.petsc",FILE_MODE_WRITE,&viewer);
405: VecView(x, viewer);
406: PetscViewerDestroy(&viewer);
407: }
409: flg = PETSC_FALSE;
410: PetscOptionsGetBool(NULL, "-ksp_reason", &flg,NULL);
411: if (flg) {
412: KSPConvergedReason reason;
413: KSPGetConvergedReason(ksp,&reason);
414: PetscPrintf(PETSC_COMM_WORLD,"KSPConvergedReason: %D\n", reason);
415: }
417: } /* while (num_numfac--) */
419: /*
420: Free work space. All PETSc objects should be destroyed when they
421: are no longer needed.
422: */
423: MatDestroy(&A); VecDestroy(&b);
424: VecDestroy(&u); VecDestroy(&x);
425: VecDestroy(&b2);
426: KSPDestroy(&ksp);
427: if (flgB) { MatDestroy(&B); }
428: PetscPreLoadEnd();
429: /* -----------------------------------------------------------
430: End of linear solver loop
431: ----------------------------------------------------------- */
433: PetscFinalize();
434: return 0;
435: }