Actual source code: ex17.c
1: static const char help[] = "Newton's method to solve a two-variable system, sequentially.\n"
2: "The same problem is solved twice - i) fully assembled system + ii) block system\n\n";
4: /*T
5: Concepts: SNES^basic uniprocessor example, block objects
6: Processors: 1
7: T*/
9: /*
10: Include "petscsnes.h" so that we can use SNES solvers. Note that this
11: file automatically includes:
12: petscsys.h - base PETSc routines petscvec.h - vectors
13: petscsys.h - system routines petscmat.h - matrices
14: petscis.h - index sets petscksp.h - Krylov subspace methods
15: petscviewer.h - viewers petscpc.h - preconditioners
16: petscksp.h - linear solvers
17: */
18: #include <petscsnes.h>
20: /*
21: This example is block version of the test found at
22: ${PETSC_DIR}/src/snes/examples/tutorials/ex1.c
23: In this test we replace the Jacobian systems
24: [J]{x} = {F}
25: where
27: [J] = ( j_00, j_01 ), {x} = ( x_0, x_1 )^T, {F} = ( f_0, f_1 )^T
28: ( j_10, j_11 )
29: where [J] \in \mathbb^{2 \times 2}, {x},{F} \in \mathbb^{2 \times 1},
31: with a block system in which each block is of length 1.
32: i.e. The block system is thus
34: [J] = ( [j00], [j01] ), {x} = ( {x0}, {x1} )^T, {F} = ( {f0}, {f1} )^T
35: ( [j10], [j11] )
36: where
37: [j00], [j01], [j10], [j11] \in \mathbb^{1 \times 1}
38: {x0}, {x1}, {f0}, {f1} \in \mathbb^{1 \times 1}
40: In practice we would not bother defing blocks of size one, and would instead assemble the
41: full system. This is just a simple test to illustrate how to manipulate the blocks and
42: to confirm the implementation is correct.
43: */
45: /*
46: User-defined routines
47: */
48: static PetscErrorCode FormJacobian1(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
49: static PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);
50: static PetscErrorCode FormJacobian2(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
51: static PetscErrorCode FormFunction2(SNES,Vec,Vec,void*);
52: static PetscErrorCode FormJacobian1_block(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
53: static PetscErrorCode FormFunction1_block(SNES,Vec,Vec,void*);
54: static PetscErrorCode FormJacobian2_block(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
55: static PetscErrorCode FormFunction2_block(SNES,Vec,Vec,void*);
60: static PetscErrorCode assembled_system(void)
61: {
62: SNES snes; /* nonlinear solver context */
63: KSP ksp; /* linear solver context */
64: PC pc; /* preconditioner context */
65: Vec x,r; /* solution, residual vectors */
66: Mat J; /* Jacobian matrix */
68: PetscInt its;
69: PetscScalar pfive = .5,*xx;
70: PetscBool flg;
73: PetscPrintf( PETSC_COMM_WORLD, "\n\n========================= Assembled system =========================\n\n" );
75: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
76: Create nonlinear solver context
77: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
79: SNESCreate(PETSC_COMM_WORLD,&snes);
81: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
82: Create matrix and vector data structures; set corresponding routines
83: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
85: /*
86: Create vectors for solution and nonlinear function
87: */
88: VecCreateSeq(PETSC_COMM_SELF,2,&x);
89: VecDuplicate(x,&r);
91: /*
92: Create Jacobian matrix data structure
93: */
94: MatCreate(PETSC_COMM_SELF,&J);
95: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
96: MatSetFromOptions(J);
98: PetscOptionsHasName(PETSC_NULL,"-hard",&flg);
99: if (!flg) {
100: /*
101: Set function evaluation routine and vector.
102: */
103: SNESSetFunction(snes,r,FormFunction1,PETSC_NULL);
105: /*
106: Set Jacobian matrix data structure and Jacobian evaluation routine
107: */
108: SNESSetJacobian(snes,J,J,FormJacobian1,PETSC_NULL);
109: } else {
110: SNESSetFunction(snes,r,FormFunction2,PETSC_NULL);
111: SNESSetJacobian(snes,J,J,FormJacobian2,PETSC_NULL);
112: }
114: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115: Customize nonlinear solver; set runtime options
116: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
118: /*
119: Set linear solver defaults for this problem. By extracting the
120: KSP, KSP, and PC contexts from the SNES context, we can then
121: directly call any KSP, KSP, and PC routines to set various options.
122: */
123: SNESGetKSP(snes,&ksp);
124: KSPGetPC(ksp,&pc);
125: PCSetType(pc,PCNONE);
126: KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);
128: /*
129: Set SNES/KSP/KSP/PC runtime options, e.g.,
130: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
131: These options will override those specified above as long as
132: SNESSetFromOptions() is called _after_ any other customization
133: routines.
134: */
135: SNESSetFromOptions(snes);
137: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138: Evaluate initial guess; then solve nonlinear system
139: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140: if (!flg) {
141: VecSet(x,pfive);
142: } else {
143: VecGetArray(x,&xx);
144: xx[0] = 2.0; xx[1] = 3.0;
145: VecRestoreArray(x,&xx);
146: }
147: /*
148: Note: The user should initialize the vector, x, with the initial guess
149: for the nonlinear solver prior to calling SNESSolve(). In particular,
150: to employ an initial guess of zero, the user should explicitly set
151: this vector to zero by calling VecSet().
152: */
154: SNESSolve(snes,PETSC_NULL,x);
155: SNESGetIterationNumber(snes,&its);
156: if (flg) {
157: Vec f;
158: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
159: SNESGetFunction(snes,&f,0,0);
160: VecView(r,PETSC_VIEWER_STDOUT_WORLD);
161: }
163: PetscPrintf(PETSC_COMM_SELF,"number of Newton iterations = %D\n\n",its);
165: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166: Free work space. All PETSc objects should be destroyed when they
167: are no longer needed.
168: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
170: VecDestroy(&x); VecDestroy(&r);
171: MatDestroy(&J); SNESDestroy(&snes);
172: return(0);
173: }
174: /* ------------------------------------------------------------------- */
177: /*
178: FormFunction1 - Evaluates nonlinear function, F(x).
180: Input Parameters:
181: . snes - the SNES context
182: . x - input vector
183: . dummy - optional user-defined context (not used here)
185: Output Parameter:
186: . f - function vector
187: */
188: static PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *dummy)
189: {
191: PetscScalar *xx,*ff;
194: /*
195: Get pointers to vector data.
196: - For default PETSc vectors, VecGetArray() returns a pointer to
197: the data array. Otherwise, the routine is implementation dependent.
198: - You MUST call VecRestoreArray() when you no longer need access to
199: the array.
200: */
201: VecGetArray(x,&xx);
202: VecGetArray(f,&ff);
204: /*
205: Compute function
206: */
207: ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
208: ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;
211: /*
212: Restore vectors
213: */
214: VecRestoreArray(x,&xx);
215: VecRestoreArray(f,&ff);
216: return(0);
217: }
218: /* ------------------------------------------------------------------- */
221: /*
222: FormJacobian1 - Evaluates Jacobian matrix.
224: Input Parameters:
225: . snes - the SNES context
226: . x - input vector
227: . dummy - optional user-defined context (not used here)
229: Output Parameters:
230: . jac - Jacobian matrix
231: . B - optionally different preconditioning matrix
232: . flag - flag indicating matrix structure
233: */
234: static PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
235: {
236: PetscScalar *xx,A[4];
238: PetscInt idx[2] = {0,1};
241: /*
242: Get pointer to vector data
243: */
244: VecGetArray(x,&xx);
246: /*
247: Compute Jacobian entries and insert into matrix.
248: - Since this is such a small problem, we set all entries for
249: the matrix at once.
250: */
251: A[0] = 2.0*xx[0] + xx[1]; A[1] = xx[0];
252: A[2] = xx[1]; A[3] = xx[0] + 2.0*xx[1];
253: MatSetValues(*jac,2,idx,2,idx,A,INSERT_VALUES);
254: *flag = SAME_NONZERO_PATTERN;
256: /*
257: Restore vector
258: */
259: VecRestoreArray(x,&xx);
261: /*
262: Assemble matrix
263: */
264: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
265: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
266: return(0);
267: }
270: /* ------------------------------------------------------------------- */
273: static PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
274: {
276: PetscScalar *xx,*ff;
279: /*
280: Get pointers to vector data.
281: - For default PETSc vectors, VecGetArray() returns a pointer to
282: the data array. Otherwise, the routine is implementation dependent.
283: - You MUST call VecRestoreArray() when you no longer need access to
284: the array.
285: */
286: VecGetArray(x,&xx);
287: VecGetArray(f,&ff);
289: /*
290: Compute function
291: */
292: ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
293: ff[1] = xx[1];
295: /*
296: Restore vectors
297: */
298: VecRestoreArray(x,&xx);
299: VecRestoreArray(f,&ff);
300: return(0);
301: }
302: /* ------------------------------------------------------------------- */
305: static PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
306: {
307: PetscScalar *xx,A[4];
309: PetscInt idx[2] = {0,1};
312: /*
313: Get pointer to vector data
314: */
315: VecGetArray(x,&xx);
317: /*
318: Compute Jacobian entries and insert into matrix.
319: - Since this is such a small problem, we set all entries for
320: the matrix at once.
321: */
322: A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
323: A[2] = 0.0; A[3] = 1.0;
324: MatSetValues(*jac,2,idx,2,idx,A,INSERT_VALUES);
325: *flag = SAME_NONZERO_PATTERN;
327: /*
328: Restore vector
329: */
330: VecRestoreArray(x,&xx);
332: /*
333: Assemble matrix
334: */
335: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
336: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
337: return(0);
338: }
342: static int block_system(void)
343: {
344: SNES snes; /* nonlinear solver context */
345: KSP ksp; /* linear solver context */
346: PC pc; /* preconditioner context */
347: Vec x,r; /* solution, residual vectors */
348: Mat J; /* Jacobian matrix */
350: PetscInt its;
351: PetscScalar pfive = .5;
352: PetscBool flg;
354: Mat j11, j12, j21, j22;
355: Vec x1, x2, r1, r2;
356: Vec bv;
357: Vec bx[2];
358: Mat bA[2][2];
361: PetscPrintf( PETSC_COMM_WORLD, "\n\n========================= Block system =========================\n\n" );
363: SNESCreate(PETSC_COMM_WORLD,&snes);
365: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
366: Create matrix and vector data structures; set corresponding routines
367: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
369: /*
370: Create sub vectors for solution and nonlinear function
371: */
372: VecCreateSeq(PETSC_COMM_SELF,1,&x1);
373: VecDuplicate(x1,&r1);
375: VecCreateSeq(PETSC_COMM_SELF,1,&x2);
376: VecDuplicate(x2,&r2);
378: /*
379: Create the block vectors
380: */
381: bx[0] = x1;
382: bx[1] = x2;
383: VecCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL,bx,&x);
384: VecAssemblyBegin(x);
385: VecAssemblyEnd(x);
386: VecDestroy(&x1);
387: VecDestroy(&x2);
389: bx[0] = r1;
390: bx[1] = r2;
391: VecCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL,bx,&r);
392: VecDestroy(&r1);
393: VecDestroy(&r2);
394: VecAssemblyBegin(r);
395: VecAssemblyEnd(r);
397: /*
398: Create sub Jacobian matrix data structure
399: */
400: MatCreate( PETSC_COMM_WORLD, &j11 );
401: MatSetSizes( j11, 1, 1, 1, 1 );
402: MatSetType( j11, MATSEQAIJ );
404: MatCreate( PETSC_COMM_WORLD, &j12 );
405: MatSetSizes( j12, 1, 1, 1, 1 );
406: MatSetType( j12, MATSEQAIJ );
408: MatCreate( PETSC_COMM_WORLD, &j21 );
409: MatSetSizes( j21, 1, 1, 1, 1 );
410: MatSetType( j21, MATSEQAIJ );
412: MatCreate( PETSC_COMM_WORLD, &j22 );
413: MatSetSizes( j22, PETSC_DECIDE, PETSC_DECIDE, 1, 1 );
414: MatSetType( j22, MATSEQAIJ );
416: /*
417: Create block Jacobian matrix data structure
418: */
419: bA[0][0] = j11;
420: bA[0][1] = j12;
421: bA[1][0] = j21;
422: bA[1][1] = j22;
423: MatCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL,2,PETSC_NULL,&bA[0][0],&J);
424: MatNestSetVecType(J,VECNEST);
425: MatDestroy(&j11);
426: MatDestroy(&j12);
427: MatDestroy(&j21);
428: MatDestroy(&j22);
430: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
431: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
433: PetscOptionsHasName(PETSC_NULL,"-hard",&flg);
434: if (!flg) {
435: /*
436: Set function evaluation routine and vector.
437: */
438: SNESSetFunction(snes,r,FormFunction1_block,PETSC_NULL);
440: /*
441: Set Jacobian matrix data structure and Jacobian evaluation routine
442: */
443: SNESSetJacobian(snes,J,J,FormJacobian1_block,PETSC_NULL);
444: } else {
445: SNESSetFunction(snes,r,FormFunction2_block,PETSC_NULL);
446: SNESSetJacobian(snes,J,J,FormJacobian2_block,PETSC_NULL);
447: }
449: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
450: Customize nonlinear solver; set runtime options
451: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
453: /*
454: Set linear solver defaults for this problem. By extracting the
455: KSP, KSP, and PC contexts from the SNES context, we can then
456: directly call any KSP, KSP, and PC routines to set various options.
457: */
458: SNESGetKSP(snes,&ksp);
459: KSPGetPC(ksp,&pc);
460: PCSetType(pc,PCNONE);
461: KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);
463: /*
464: Set SNES/KSP/KSP/PC runtime options, e.g.,
465: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
466: These options will override those specified above as long as
467: SNESSetFromOptions() is called _after_ any other customization
468: routines.
469: */
470: SNESSetFromOptions(snes);
472: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
473: Evaluate initial guess; then solve nonlinear system
474: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
475: if (!flg) {
476: VecSet(x,pfive);
477: } else {
478: Vec *vecs;
479: VecNestGetSubVecs( x, PETSC_NULL, &vecs );
480: bv = vecs[0];
481: // VecBlockGetSubVec( x, 0, &bv );
482: VecSetValue( bv, 0, 2.0, INSERT_VALUES ); /* xx[0] = 2.0; */
483: VecAssemblyBegin(bv);
484: VecAssemblyEnd(bv);
486: // VecBlockGetSubVec( x, 1, &bv );
487: bv = vecs[1];
488: VecSetValue( bv, 0, 3.0, INSERT_VALUES ); /* xx[1] = 3.0; */
489: VecAssemblyBegin(bv);
490: VecAssemblyEnd(bv);
491: }
492: /*
493: Note: The user should initialize the vector, x, with the initial guess
494: for the nonlinear solver prior to calling SNESSolve(). In particular,
495: to employ an initial guess of zero, the user should explicitly set
496: this vector to zero by calling VecSet().
497: */
498: SNESSolve(snes,PETSC_NULL,x);
499: SNESGetIterationNumber(snes,&its);
500: if (flg) {
501: Vec f;
502: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
503: SNESGetFunction(snes,&f,0,0);
504: VecView(r,PETSC_VIEWER_STDOUT_WORLD);
505: }
507: PetscPrintf(PETSC_COMM_SELF,"number of Newton iterations = %D\n\n",its);
509: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
510: Free work space. All PETSc objects should be destroyed when they
511: are no longer needed.
512: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
513: VecDestroy(&x); VecDestroy(&r);
514: MatDestroy(&J); SNESDestroy(&snes);
516: return(0);
517: }
518: /* ------------------------------------------------------------------- */
521: static PetscErrorCode FormFunction1_block(SNES snes,Vec x,Vec f,void *dummy)
522: {
523: Vec *xx, *ff, x1,x2, f1,f2;
524: PetscScalar ff_0, ff_1;
525: PetscScalar xx_0, xx_1;
526: PetscInt index,nb;
530: /* get blocks for function */
531: VecNestGetSubVecs( f, &nb, &ff );
532: f1 = ff[0]; f2 = ff[1];
534: /* get blocks for solution */
535: VecNestGetSubVecs( x, &nb, &xx );
536: x1 = xx[0]; x2 = xx[1];
538: /* get solution values */
539: index = 0;
540: VecGetValues( x1,1, &index, &xx_0 );
541: VecGetValues( x2,1, &index, &xx_1 );
543: /* Compute function */
544: ff_0 = xx_0*xx_0 + xx_0*xx_1 - 3.0;
545: ff_1 = xx_0*xx_1 + xx_1*xx_1 - 6.0;
547: /* set function values */
548: VecSetValue( f1, index, ff_0, INSERT_VALUES );
550: VecSetValue( f2, index, ff_1, INSERT_VALUES );
552: VecAssemblyBegin(f);
553: VecAssemblyEnd(f);
555: return(0);
556: }
557: /* ------------------------------------------------------------------- */
560: static PetscErrorCode FormJacobian1_block(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
561: {
562: Vec *xx, x1,x2;
563: PetscScalar xx_0, xx_1;
564: PetscInt index,nb;
565: PetscScalar A_00, A_01, A_10, A_11;
566: Mat j11, j12, j21, j22;
567: Mat **mats;
571: /* get blocks for solution */
572: VecNestGetSubVecs( x, &nb, &xx );
573: x1 = xx[0]; x2 = xx[1];
575: /* get solution values */
576: index = 0;
577: VecGetValues( x1,1, &index, &xx_0 );
578: VecGetValues( x2,1, &index, &xx_1 );
580: /* get block matrices */
581: MatNestGetSubMats(*jac,PETSC_NULL,PETSC_NULL,&mats);
582: j11 = mats[0][0];
583: j12 = mats[0][1];
584: j21 = mats[1][0];
585: j22 = mats[1][1];
587: /* compute jacobian entries */
588: A_00 = 2.0*xx_0 + xx_1;
589: A_01 = xx_0;
590: A_10 = xx_1;
591: A_11 = xx_0 + 2.0*xx_1;
593: /* set jacobian values */
594: MatSetValue( j11, 0,0, A_00, INSERT_VALUES);
595: MatSetValue( j12, 0,0, A_01, INSERT_VALUES);
596: MatSetValue( j21, 0,0, A_10, INSERT_VALUES);
597: MatSetValue( j22, 0,0, A_11, INSERT_VALUES);
599: *flag = SAME_NONZERO_PATTERN;
601: /* Assemble sub matrix */
602: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
603: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
604: return(0);
605: }
608: /* ------------------------------------------------------------------- */
611: static PetscErrorCode FormFunction2_block(SNES snes,Vec x,Vec f,void *dummy)
612: {
614: PetscScalar *xx,*ff;
617: /*
618: Get pointers to vector data.
619: - For default PETSc vectors, VecGetArray() returns a pointer to
620: the data array. Otherwise, the routine is implementation dependent.
621: - You MUST call VecRestoreArray() when you no longer need access to
622: the array.
623: */
624: VecGetArray(x,&xx);
625: VecGetArray(f,&ff);
627: /*
628: Compute function
629: */
630: ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
631: ff[1] = xx[1];
633: /*
634: Restore vectors
635: */
636: VecRestoreArray(x,&xx);
637: VecRestoreArray(f,&ff);
638: return(0);
639: }
640: /* ------------------------------------------------------------------- */
643: static PetscErrorCode FormJacobian2_block(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
644: {
645: PetscScalar *xx,A[4];
647: PetscInt idx[2] = {0,1};
650: /*
651: Get pointer to vector data
652: */
653: VecGetArray(x,&xx);
655: /*
656: Compute Jacobian entries and insert into matrix.
657: - Since this is such a small problem, we set all entries for
658: the matrix at once.
659: */
660: A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
661: A[2] = 0.0; A[3] = 1.0;
662: MatSetValues(*jac,2,idx,2,idx,A,INSERT_VALUES);
663: *flag = SAME_NONZERO_PATTERN;
665: /*
666: Restore vector
667: */
668: VecRestoreArray(x,&xx);
670: /*
671: Assemble matrix
672: */
673: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
674: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
675: return(0);
676: }
680: int main(int argc,char **argv)
681: {
682: PetscMPIInt size;
685: PetscInitialize(&argc,&argv,(char *)0,help);
687: MPI_Comm_size(PETSC_COMM_WORLD,&size);
688: if (size != 1) SETERRQ(PETSC_COMM_WORLD, 1,"This is a uniprocessor example only!");
690: assembled_system();
692: block_system();
694: PetscFinalize();
695: return 0;
696: }