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