1 
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*/
27 
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>
38 
39 #undef __FUNCT__
40 #define __FUNCT__ "main"
main(int argc,char ** args)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;
50   PetscErrorCode ierr;
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];
57 
58   PetscInitialize(&argc,&args,(char*)0,help);
59   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
60   ierr = PetscOptionsGetBool(NULL,"-table",&table,NULL);CHKERRQ(ierr);
61   ierr = PetscOptionsGetBool(NULL,"-trans",&trans,NULL);CHKERRQ(ierr);
62   ierr = PetscOptionsGetBool(NULL,"-initialguess",&initialguess,NULL);CHKERRQ(ierr);
63   ierr = PetscOptionsGetBool(NULL,"-output_solution",&outputSoln,NULL);CHKERRQ(ierr);
64   ierr = PetscOptionsGetString(NULL,"-initialguessfilename",initialguessfilename,PETSC_MAX_PATH_LEN,&initialguessfile);CHKERRQ(ierr);
65 
66   /*
67      Determine files from which we read the two linear systems
68      (matrix and right-hand-side vector).
69   */
70   ierr = PetscOptionsGetString(NULL,"-f",file[0],PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
71   if (flg) {
72     ierr    = PetscStrcpy(file[1],file[0]);CHKERRQ(ierr);
73     preload = PETSC_FALSE;
74   } else {
75     ierr = PetscOptionsGetString(NULL,"-f0",file[0],PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
76     if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f0 or -f option");
77     ierr = PetscOptionsGetString(NULL,"-f1",file[1],PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
78     if (!flg) preload = PETSC_FALSE;   /* don't bother with second system */
79   }
80 
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");
94 
95   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
96                          Load system
97    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
98 
99   /*
100      Open binary file.  Note that we use FILE_MODE_READ to indicate
101      reading from this file.
102   */
103   ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[PetscPreLoadIt],FILE_MODE_READ,&fd);CHKERRQ(ierr);
104 
105   /*
106      Load the matrix and vector; then destroy the viewer.
107   */
108   ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
109   ierr = MatSetFromOptions(A);CHKERRQ(ierr);
110   ierr = MatLoad(A,fd);CHKERRQ(ierr);
111 
112   flg  = PETSC_FALSE;
113   ierr = PetscOptionsGetString(NULL,"-rhs",file[2],PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
114   ierr = VecCreate(PETSC_COMM_WORLD,&b);CHKERRQ(ierr);
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       ierr = PetscInfo(0,"Using vector of ones for RHS\n");CHKERRQ(ierr);
120       ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr);
121       ierr = VecSetSizes(b,m,PETSC_DECIDE);CHKERRQ(ierr);
122       ierr = VecSetFromOptions(b);CHKERRQ(ierr);
123       ierr = VecSet(b,one);CHKERRQ(ierr);
124     } else {
125       ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);
126       ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2],FILE_MODE_READ,&fd);CHKERRQ(ierr);
127       ierr = VecSetFromOptions(b);CHKERRQ(ierr);
128       ierr = VecLoad(b,fd);CHKERRQ(ierr);
129     }
130   } else {   /* rhs is stored in the same file as matrix */
131     ierr = VecSetFromOptions(b);CHKERRQ(ierr);
132     ierr = VecLoad(b,fd);CHKERRQ(ierr);
133   }
134   ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);
135 
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   ierr = PetscOptionsGetBool(NULL, "-test_zeropivot", &flg,NULL);CHKERRQ(ierr);
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     ierr = MatGetRow(A,row,&ncols,&cols,&vals);CHKERRQ(ierr);
148     ierr = PetscMalloc(sizeof(PetscScalar)*(ncols+1),&zeros);
149     ierr = PetscMemzero(zeros,(ncols+1)*sizeof(PetscScalar));CHKERRQ(ierr);
150     ierr = PetscOptionsGetBool(NULL, "-set_row_zero", &flg1,NULL);CHKERRQ(ierr);
151     if (flg1) {   /* set entire row as zero */
152       ierr = MatSetValues(A,1,&row,ncols,cols,zeros,INSERT_VALUES);CHKERRQ(ierr);
153     } else {   /* only set (row,row) entry as zero */
154       ierr = MatSetValues(A,1,&row,1,&row,zeros,INSERT_VALUES);CHKERRQ(ierr);
155     }
156     ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
157     ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
158   }
159 
160   /* Check whether A is symmetric */
161   flg  = PETSC_FALSE;
162   ierr = PetscOptionsGetBool(NULL, "-check_symmetry", &flg,NULL);CHKERRQ(ierr);
163   if (flg) {
164     Mat Atrans;
165     ierr = MatTranspose(A, MAT_INITIAL_MATRIX,&Atrans);
166     ierr = MatEqual(A, Atrans, &isSymmetric);
167     if (isSymmetric) {
168       ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
169     } else {
170       PetscPrintf(PETSC_COMM_WORLD,"Warning: A is non-symmetric \n");CHKERRQ(ierr);
171     }
172     ierr = MatDestroy(&Atrans);CHKERRQ(ierr);
173   }
174 
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   */
179 
180   ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr);
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   ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr);
183   ierr = VecGetSize(b,&m);CHKERRQ(ierr);
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;
188 
189     ierr = VecCreate(PETSC_COMM_WORLD,&tmp);CHKERRQ(ierr);
190     ierr = VecSetSizes(tmp,n,PETSC_DECIDE);CHKERRQ(ierr);
191     ierr = VecSetFromOptions(tmp);CHKERRQ(ierr);
192     ierr = VecGetOwnershipRange(b,&start,&end);CHKERRQ(ierr);
193     ierr = VecGetLocalSize(b,&mvec);CHKERRQ(ierr);
194     ierr = VecGetArray(b,&bold);CHKERRQ(ierr);
195     for (j=0; j<mvec; j++) {
196       indx = start+j;
197       ierr = VecSetValues(tmp,1,&indx,bold+j,INSERT_VALUES);CHKERRQ(ierr);
198     }
199     ierr = VecRestoreArray(b,&bold);CHKERRQ(ierr);
200     ierr = VecDestroy(&b);CHKERRQ(ierr);
201     ierr = VecAssemblyBegin(tmp);CHKERRQ(ierr);
202     ierr = VecAssemblyEnd(tmp);CHKERRQ(ierr);
203     b    = tmp;
204   }
205 
206   ierr = MatGetVecs(A,&x,NULL);CHKERRQ(ierr);
207   ierr = VecDuplicate(b,&u);CHKERRQ(ierr);
208   if (initialguessfile) {
209     PetscViewer viewer2;
210     ierr         = PetscViewerBinaryOpen(PETSC_COMM_WORLD,initialguessfilename,FILE_MODE_READ,&viewer2);CHKERRQ(ierr);
211     ierr         = VecLoad(x,viewer2);CHKERRQ(ierr);
212     ierr         = PetscViewerDestroy(&viewer2);CHKERRQ(ierr);
213     initialguess = PETSC_TRUE;
214   } else if (initialguess) {
215     ierr = VecSet(x,1.0);CHKERRQ(ierr);
216   } else {
217     ierr = VecSet(x,0.0);CHKERRQ(ierr);
218   }
219 
220 
221   /* Check scaling in A */
222   flg  = PETSC_FALSE;
223   ierr = PetscOptionsGetBool(NULL, "-check_scaling", &flg,NULL);CHKERRQ(ierr);
224   if (flg) {
225     Vec       max, min;
226     PetscInt  idx;
227     PetscReal val;
228 
229     ierr = VecDuplicate(x, &max);CHKERRQ(ierr);
230     ierr = VecDuplicate(x, &min);CHKERRQ(ierr);
231     ierr = MatGetRowMaxAbs(A, max, NULL);CHKERRQ(ierr);
232     ierr = MatGetRowMinAbs(A, min, NULL);CHKERRQ(ierr);
233     {
234       PetscViewer viewer;
235 
236       ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, "max.data", &viewer);CHKERRQ(ierr);
237       ierr = VecView(max, viewer);CHKERRQ(ierr);
238       ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
239       ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, "min.data", &viewer);CHKERRQ(ierr);
240       ierr = VecView(min, viewer);CHKERRQ(ierr);
241       ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
242     }
243     ierr = VecView(max, PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);
244     ierr = VecMax(max, &idx, &val);CHKERRQ(ierr);
245     ierr = PetscPrintf(PETSC_COMM_WORLD, "Largest max row element %G at row %d\n", val, idx);CHKERRQ(ierr);
246     ierr = VecView(min, PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);
247     ierr = VecMin(min, &idx, &val);CHKERRQ(ierr);
248     ierr = PetscPrintf(PETSC_COMM_WORLD, "Smallest min row element %G at row %d\n", val, idx);CHKERRQ(ierr);
249     ierr = VecMin(max, &idx, &val);CHKERRQ(ierr);
250     ierr = PetscPrintf(PETSC_COMM_WORLD, "Smallest max row element %G at row %d\n", val, idx);CHKERRQ(ierr);
251     ierr = VecPointwiseDivide(max, max, min);CHKERRQ(ierr);
252     ierr = VecMax(max, &idx, &val);CHKERRQ(ierr);
253     ierr = PetscPrintf(PETSC_COMM_WORLD, "Largest row ratio %G at row %d\n", val, idx);CHKERRQ(ierr);
254     ierr = VecView(max, PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);
255     ierr = VecDestroy(&max);CHKERRQ(ierr);
256     ierr = VecDestroy(&min);CHKERRQ(ierr);
257   }
258 
259   /*  ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */
260   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
261                     Setup solve for system
262    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
263   /*
264      Conclude profiling last stage; begin profiling next stage.
265   */
266   PetscPreLoadStage("KSPSetUpSolve");
267 
268   /*
269      We also explicitly time this stage via PetscTime()
270   */
271   ierr = PetscTime(&tsetup1);CHKERRQ(ierr);
272 
273   /*
274      Create linear solver; set operators; set runtime options.
275   */
276   ierr       = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
277   ierr       = KSPSetInitialGuessNonzero(ksp,initialguess);CHKERRQ(ierr);
278   num_numfac = 1;
279   ierr       = PetscOptionsGetInt(NULL,"-num_numfac",&num_numfac,NULL);CHKERRQ(ierr);
280   while (num_numfac--) {
281     PetscBool lsqr;
282     char      str[32];
283     ierr = PetscOptionsGetString(NULL,"-ksp_type",str,32,&lsqr);CHKERRQ(ierr);
284     if (lsqr) {
285       ierr = PetscStrcmp("lsqr",str,&lsqr);CHKERRQ(ierr);
286     }
287     if (lsqr) {
288       Mat BtB;
289       ierr = MatTransposeMatMult(A,A,MAT_INITIAL_MATRIX,4,&BtB);CHKERRQ(ierr);
290       ierr = KSPSetOperators(ksp,A,BtB,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
291       ierr = MatDestroy(&BtB);CHKERRQ(ierr);
292     } else {
293       ierr = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
294     }
295     ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
296 
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     ierr   = KSPSetUp(ksp);CHKERRQ(ierr);
304     ierr   = KSPSetUpOnBlocks(ksp);CHKERRQ(ierr);
305     ierr   = PetscTime(&tsetup2);CHKERRQ(ierr);
306     tsetup = tsetup2 - tsetup1;
307 
308     /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
309                          Solve system
310       - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
311 
312     /*
313      Solve linear system; we also explicitly time this stage.
314     */
315     ierr = PetscTime(&tsolve1);CHKERRQ(ierr);
316     if (trans) {
317       ierr = KSPSolveTranspose(ksp,b,x);CHKERRQ(ierr);
318       ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
319     } else {
320       PetscInt num_rhs=1;
321       ierr   = PetscOptionsGetInt(NULL,"-num_rhs",&num_rhs,NULL);CHKERRQ(ierr);
322       cknorm = PETSC_FALSE;
323       ierr   = PetscOptionsGetBool(NULL,"-cknorm",&cknorm,NULL);CHKERRQ(ierr);
324       while (num_rhs--) {
325         if (num_rhs == 1) VecSet(x,0.0);
326         ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
327       }
328       ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
329       if (cknorm) {     /* Check error for each rhs */
330         if (trans) {
331           ierr = MatMultTranspose(A,x,u);CHKERRQ(ierr);
332         } else {
333           ierr = MatMult(A,x,u);CHKERRQ(ierr);
334         }
335         ierr = VecAXPY(u,-1.0,b);CHKERRQ(ierr);
336         ierr = VecNorm(u,NORM_2,&norm);CHKERRQ(ierr);
337         ierr = PetscPrintf(PETSC_COMM_WORLD,"  Number of iterations = %3D\n",its);CHKERRQ(ierr);
338         if (norm < 1.e-12) {
339           ierr = PetscPrintf(PETSC_COMM_WORLD,"  Residual norm < 1.e-12\n");CHKERRQ(ierr);
340         } else {
341           ierr = PetscPrintf(PETSC_COMM_WORLD,"  Residual norm %G\n",norm);CHKERRQ(ierr);
342         }
343       }
344     }   /* while (num_rhs--) */
345     ierr   = PetscTime(&tsolve2);CHKERRQ(ierr);
346     tsolve = tsolve2 - tsolve1;
347 
348     /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
349           Check error, print output, free data structures.
350      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
351 
352     /*
353        Check error
354     */
355     if (trans) {
356       ierr = MatMultTranspose(A,x,u);CHKERRQ(ierr);
357     } else {
358       ierr = MatMult(A,x,u);CHKERRQ(ierr);
359     }
360     ierr = VecAXPY(u,-1.0,b);CHKERRQ(ierr);
361     ierr = VecNorm(u,NORM_2,&norm);CHKERRQ(ierr);
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;
371 
372       /*
373        Open a string viewer; then write info to it.
374       */
375       ierr = PetscViewerStringOpen(PETSC_COMM_WORLD,kspinfo,120,&viewer);CHKERRQ(ierr);
376       ierr = KSPView(ksp,viewer);CHKERRQ(ierr);
377       ierr = PetscStrrchr(file[PetscPreLoadIt],'/',&matrixname);CHKERRQ(ierr);
378       ierr = 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);CHKERRQ(ierr);
380 
381       /*
382         Destroy the viewer
383       */
384       ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
385     } else {
386       ierr = PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);CHKERRQ(ierr);
387       if (norm < 1.e-12) {
388         ierr = PetscPrintf(PETSC_COMM_WORLD,"  Residual norm < 1.e-12\n");CHKERRQ(ierr);
389       } else {
390         ierr = PetscPrintf(PETSC_COMM_WORLD,"Residual norm %G\n",norm);CHKERRQ(ierr);
391       }
392     }
393     ierr = PetscOptionsGetString(NULL,"-solution",file[3],PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
394     if (flg) {
395       PetscViewer viewer;
396       Vec         xstar;
397       PetscReal   norm;
398 
399       ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[3],FILE_MODE_READ,&viewer);CHKERRQ(ierr);
400       ierr = VecCreate(PETSC_COMM_WORLD,&xstar);CHKERRQ(ierr);
401       ierr = VecLoad(xstar,viewer);CHKERRQ(ierr);
402       ierr = VecAXPY(xstar, -1.0, x);CHKERRQ(ierr);
403       ierr = VecNorm(xstar, NORM_2, &norm);CHKERRQ(ierr);
404       ierr = PetscPrintf(PETSC_COMM_WORLD, "Error norm %G\n", norm);CHKERRQ(ierr);
405       ierr = VecDestroy(&xstar);CHKERRQ(ierr);
406       ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
407     }
408     if (outputSoln) {
409       PetscViewer viewer;
410 
411       ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"solution.petsc",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
412       ierr = VecView(x, viewer);CHKERRQ(ierr);
413       ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
414     }
415 
416     flg  = PETSC_FALSE;
417     ierr = PetscOptionsGetBool(NULL, "-ksp_reason", &flg,NULL);CHKERRQ(ierr);
418     if (flg) {
419       KSPConvergedReason reason;
420       ierr = KSPGetConvergedReason(ksp,&reason);CHKERRQ(ierr);
421       PetscPrintf(PETSC_COMM_WORLD,"KSPConvergedReason: %D\n", reason);
422     }
423 
424   }   /* while (num_numfac--) */
425 
426   /*
427      Free work space.  All PETSc objects should be destroyed when they
428      are no longer needed.
429   */
430   ierr = MatDestroy(&A);CHKERRQ(ierr); ierr = VecDestroy(&b);CHKERRQ(ierr);
431   ierr = VecDestroy(&u);CHKERRQ(ierr); ierr = VecDestroy(&x);CHKERRQ(ierr);
432   ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
433   PetscPreLoadEnd();
434   /* -----------------------------------------------------------
435                       End of linear solver loop
436      ----------------------------------------------------------- */
437 
438   ierr = PetscFinalize();
439   return 0;
440 }
441