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: /*
5: Include "petscsnes.h" so that we can use SNES solvers. Note that this
6: file automatically includes:
7: petscsys.h - base PETSc routines petscvec.h - vectors
8: petscsys.h - system routines petscmat.h - matrices
9: petscis.h - index sets petscksp.h - Krylov subspace methods
10: petscviewer.h - viewers petscpc.h - preconditioners
11: petscksp.h - linear solvers
12: */
13: #include <petscsnes.h>
15: /*
16: This example is block version of the test found at
17: ${PETSC_DIR}/src/snes/tutorials/ex1.c
18: In this test we replace the Jacobian systems
19: [J]{x} = {F}
20: where
22: [J] = (j_00, j_01), {x} = (x_0, x_1)^T, {F} = (f_0, f_1)^T
23: (j_10, j_11)
24: where [J] \in \mathbb^{2 \times 2}, {x},{F} \in \mathbb^{2 \times 1},
26: with a block system in which each block is of length 1.
27: i.e. The block system is thus
29: [J] = ([j00], [j01]), {x} = ({x0}, {x1})^T, {F} = ({f0}, {f1})^T
30: ([j10], [j11])
31: where
32: [j00], [j01], [j10], [j11] \in \mathbb^{1 \times 1}
33: {x0}, {x1}, {f0}, {f1} \in \mathbb^{1 \times 1}
35: In practice we would not bother defing blocks of size one, and would instead assemble the
36: full system. This is just a simple test to illustrate how to manipulate the blocks and
37: to confirm the implementation is correct.
38: */
40: /*
41: User-defined routines
42: */
43: static PetscErrorCode FormJacobian1(SNES, Vec, Mat, Mat, void *);
44: static PetscErrorCode FormFunction1(SNES, Vec, Vec, void *);
45: static PetscErrorCode FormJacobian2(SNES, Vec, Mat, Mat, void *);
46: static PetscErrorCode FormFunction2(SNES, Vec, Vec, void *);
47: static PetscErrorCode FormJacobian1_block(SNES, Vec, Mat, Mat, void *);
48: static PetscErrorCode FormFunction1_block(SNES, Vec, Vec, void *);
49: static PetscErrorCode FormJacobian2_block(SNES, Vec, Mat, Mat, void *);
50: static PetscErrorCode FormFunction2_block(SNES, Vec, Vec, void *);
52: static PetscErrorCode assembled_system(void)
53: {
54: SNES snes; /* nonlinear solver context */
55: KSP ksp; /* linear solver context */
56: PC pc; /* preconditioner context */
57: Vec x, r; /* solution, residual vectors */
58: Mat J; /* Jacobian matrix */
59: PetscInt its;
60: PetscScalar pfive = .5, *xx;
61: PetscBool flg;
63: PetscFunctionBeginUser;
64: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n\n========================= Assembled system =========================\n\n"));
66: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
67: Create nonlinear solver context
68: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
70: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
72: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
73: Create matrix and vector data structures; set corresponding routines
74: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
76: /*
77: Create vectors for solution and nonlinear function
78: */
79: PetscCall(VecCreateSeq(PETSC_COMM_SELF, 2, &x));
80: PetscCall(VecDuplicate(x, &r));
82: /*
83: Create Jacobian matrix data structure
84: */
85: PetscCall(MatCreate(PETSC_COMM_SELF, &J));
86: PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, 2, 2));
87: PetscCall(MatSetFromOptions(J));
88: PetscCall(MatSetUp(J));
90: PetscCall(PetscOptionsHasName(NULL, NULL, "-hard", &flg));
91: if (!flg) {
92: /*
93: Set function evaluation routine and vector.
94: */
95: PetscCall(SNESSetFunction(snes, r, FormFunction1, NULL));
97: /*
98: Set Jacobian matrix data structure and Jacobian evaluation routine
99: */
100: PetscCall(SNESSetJacobian(snes, J, J, FormJacobian1, NULL));
101: } else {
102: PetscCall(SNESSetFunction(snes, r, FormFunction2, NULL));
103: PetscCall(SNESSetJacobian(snes, J, J, FormJacobian2, NULL));
104: }
106: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: Customize nonlinear solver; set runtime options
108: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110: /*
111: Set linear solver defaults for this problem. By extracting the
112: KSP, KSP, and PC contexts from the SNES context, we can then
113: directly call any KSP, KSP, and PC routines to set various options.
114: */
115: PetscCall(SNESGetKSP(snes, &ksp));
116: PetscCall(KSPGetPC(ksp, &pc));
117: PetscCall(PCSetType(pc, PCNONE));
118: PetscCall(KSPSetTolerances(ksp, 1.e-4, PETSC_DEFAULT, PETSC_DEFAULT, 20));
120: /*
121: Set SNES/KSP/KSP/PC runtime options, e.g.,
122: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
123: These options will override those specified above as long as
124: SNESSetFromOptions() is called _after_ any other customization
125: routines.
126: */
127: PetscCall(SNESSetFromOptions(snes));
129: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130: Evaluate initial guess; then solve nonlinear system
131: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
132: if (!flg) {
133: PetscCall(VecSet(x, pfive));
134: } else {
135: PetscCall(VecGetArray(x, &xx));
136: xx[0] = 2.0;
137: xx[1] = 3.0;
138: PetscCall(VecRestoreArray(x, &xx));
139: }
140: /*
141: Note: The user should initialize the vector, x, with the initial guess
142: for the nonlinear solver prior to calling SNESSolve(). In particular,
143: to employ an initial guess of zero, the user should explicitly set
144: this vector to zero by calling VecSet().
145: */
147: PetscCall(SNESSolve(snes, NULL, x));
148: PetscCall(SNESGetIterationNumber(snes, &its));
149: if (flg) {
150: Vec f;
151: PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
152: PetscCall(SNESGetFunction(snes, &f, 0, 0));
153: PetscCall(VecView(r, PETSC_VIEWER_STDOUT_WORLD));
154: }
155: PetscCall(PetscPrintf(PETSC_COMM_SELF, "number of SNES iterations = %" PetscInt_FMT "\n\n", its));
157: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158: Free work space. All PETSc objects should be destroyed when they
159: are no longer needed.
160: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
162: PetscCall(VecDestroy(&x));
163: PetscCall(VecDestroy(&r));
164: PetscCall(MatDestroy(&J));
165: PetscCall(SNESDestroy(&snes));
166: PetscFunctionReturn(PETSC_SUCCESS);
167: }
169: /*
170: FormFunction1 - Evaluates nonlinear function, F(x).
172: Input Parameters:
173: . snes - the SNES context
174: . x - input vector
175: . dummy - optional user-defined context (not used here)
177: Output Parameter:
178: . f - function vector
179: */
180: static PetscErrorCode FormFunction1(SNES snes, Vec x, Vec f, void *dummy)
181: {
182: const PetscScalar *xx;
183: PetscScalar *ff;
185: PetscFunctionBeginUser;
186: /*
187: Get pointers to vector data.
188: - For default PETSc vectors, VecGetArray() returns a pointer to
189: the data array. Otherwise, the routine is implementation dependent.
190: - You MUST call VecRestoreArray() when you no longer need access to
191: the array.
192: */
193: PetscCall(VecGetArrayRead(x, &xx));
194: PetscCall(VecGetArray(f, &ff));
196: /*
197: Compute function
198: */
199: ff[0] = xx[0] * xx[0] + xx[0] * xx[1] - 3.0;
200: ff[1] = xx[0] * xx[1] + xx[1] * xx[1] - 6.0;
202: /*
203: Restore vectors
204: */
205: PetscCall(VecRestoreArrayRead(x, &xx));
206: PetscCall(VecRestoreArray(f, &ff));
207: PetscFunctionReturn(PETSC_SUCCESS);
208: }
210: /*
211: FormJacobian1 - Evaluates Jacobian matrix.
213: Input Parameters:
214: . snes - the SNES context
215: . x - input vector
216: . dummy - optional user-defined context (not used here)
218: Output Parameters:
219: . jac - Jacobian matrix
220: . B - optionally different preconditioning matrix
221: . flag - flag indicating matrix structure
222: */
223: static PetscErrorCode FormJacobian1(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
224: {
225: const PetscScalar *xx;
226: PetscScalar A[4];
227: PetscInt idx[2] = {0, 1};
229: PetscFunctionBeginUser;
230: /*
231: Get pointer to vector data
232: */
233: PetscCall(VecGetArrayRead(x, &xx));
235: /*
236: Compute Jacobian entries and insert into matrix.
237: - Since this is such a small problem, we set all entries for
238: the matrix at once.
239: */
240: A[0] = 2.0 * xx[0] + xx[1];
241: A[1] = xx[0];
242: A[2] = xx[1];
243: A[3] = xx[0] + 2.0 * xx[1];
244: PetscCall(MatSetValues(jac, 2, idx, 2, idx, A, INSERT_VALUES));
246: /*
247: Restore vector
248: */
249: PetscCall(VecRestoreArrayRead(x, &xx));
251: /*
252: Assemble matrix
253: */
254: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
255: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
256: PetscFunctionReturn(PETSC_SUCCESS);
257: }
259: static PetscErrorCode FormFunction2(SNES snes, Vec x, Vec f, void *dummy)
260: {
261: const PetscScalar *xx;
262: PetscScalar *ff;
264: PetscFunctionBeginUser;
265: /*
266: Get pointers to vector data.
267: - For default PETSc vectors, VecGetArray() returns a pointer to
268: the data array. Otherwise, the routine is implementation dependent.
269: - You MUST call VecRestoreArray() when you no longer need access to
270: the array.
271: */
272: PetscCall(VecGetArrayRead(x, &xx));
273: PetscCall(VecGetArray(f, &ff));
275: /*
276: Compute function
277: */
278: ff[0] = PetscSinScalar(3.0 * xx[0]) + xx[0];
279: ff[1] = xx[1];
281: /*
282: Restore vectors
283: */
284: PetscCall(VecRestoreArrayRead(x, &xx));
285: PetscCall(VecRestoreArray(f, &ff));
286: PetscFunctionReturn(PETSC_SUCCESS);
287: }
289: static PetscErrorCode FormJacobian2(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
290: {
291: const PetscScalar *xx;
292: PetscScalar A[4];
293: PetscInt idx[2] = {0, 1};
295: PetscFunctionBeginUser;
296: /*
297: Get pointer to vector data
298: */
299: PetscCall(VecGetArrayRead(x, &xx));
301: /*
302: Compute Jacobian entries and insert into matrix.
303: - Since this is such a small problem, we set all entries for
304: the matrix at once.
305: */
306: A[0] = 3.0 * PetscCosScalar(3.0 * xx[0]) + 1.0;
307: A[1] = 0.0;
308: A[2] = 0.0;
309: A[3] = 1.0;
310: PetscCall(MatSetValues(jac, 2, idx, 2, idx, A, INSERT_VALUES));
312: /*
313: Restore vector
314: */
315: PetscCall(VecRestoreArrayRead(x, &xx));
317: /*
318: Assemble matrix
319: */
320: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
321: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
322: PetscFunctionReturn(PETSC_SUCCESS);
323: }
325: static PetscErrorCode block_system(void)
326: {
327: SNES snes; /* nonlinear solver context */
328: KSP ksp; /* linear solver context */
329: PC pc; /* preconditioner context */
330: Vec x, r; /* solution, residual vectors */
331: Mat J; /* Jacobian matrix */
332: PetscInt its;
333: PetscScalar pfive = .5;
334: PetscBool flg;
336: Mat j11, j12, j21, j22;
337: Vec x1, x2, r1, r2;
338: Vec bv;
339: Vec bx[2];
340: Mat bA[2][2];
342: PetscFunctionBeginUser;
343: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n\n========================= Block system =========================\n\n"));
345: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
347: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
348: Create matrix and vector data structures; set corresponding routines
349: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
351: /*
352: Create sub vectors for solution and nonlinear function
353: */
354: PetscCall(VecCreateSeq(PETSC_COMM_SELF, 1, &x1));
355: PetscCall(VecDuplicate(x1, &r1));
357: PetscCall(VecCreateSeq(PETSC_COMM_SELF, 1, &x2));
358: PetscCall(VecDuplicate(x2, &r2));
360: /*
361: Create the block vectors
362: */
363: bx[0] = x1;
364: bx[1] = x2;
365: PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, NULL, bx, &x));
366: PetscCall(VecAssemblyBegin(x));
367: PetscCall(VecAssemblyEnd(x));
368: PetscCall(VecDestroy(&x1));
369: PetscCall(VecDestroy(&x2));
371: bx[0] = r1;
372: bx[1] = r2;
373: PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, NULL, bx, &r));
374: PetscCall(VecDestroy(&r1));
375: PetscCall(VecDestroy(&r2));
376: PetscCall(VecAssemblyBegin(r));
377: PetscCall(VecAssemblyEnd(r));
379: /*
380: Create sub Jacobian matrix data structure
381: */
382: PetscCall(MatCreate(PETSC_COMM_WORLD, &j11));
383: PetscCall(MatSetSizes(j11, 1, 1, 1, 1));
384: PetscCall(MatSetType(j11, MATSEQAIJ));
385: PetscCall(MatSetUp(j11));
387: PetscCall(MatCreate(PETSC_COMM_WORLD, &j12));
388: PetscCall(MatSetSizes(j12, 1, 1, 1, 1));
389: PetscCall(MatSetType(j12, MATSEQAIJ));
390: PetscCall(MatSetUp(j12));
392: PetscCall(MatCreate(PETSC_COMM_WORLD, &j21));
393: PetscCall(MatSetSizes(j21, 1, 1, 1, 1));
394: PetscCall(MatSetType(j21, MATSEQAIJ));
395: PetscCall(MatSetUp(j21));
397: PetscCall(MatCreate(PETSC_COMM_WORLD, &j22));
398: PetscCall(MatSetSizes(j22, PETSC_DECIDE, PETSC_DECIDE, 1, 1));
399: PetscCall(MatSetType(j22, MATSEQAIJ));
400: PetscCall(MatSetUp(j22));
401: /*
402: Create block Jacobian matrix data structure
403: */
404: bA[0][0] = j11;
405: bA[0][1] = j12;
406: bA[1][0] = j21;
407: bA[1][1] = j22;
409: PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, &bA[0][0], &J));
410: PetscCall(MatSetUp(J));
411: PetscCall(MatNestSetVecType(J, VECNEST));
412: PetscCall(MatDestroy(&j11));
413: PetscCall(MatDestroy(&j12));
414: PetscCall(MatDestroy(&j21));
415: PetscCall(MatDestroy(&j22));
417: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
418: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
420: PetscCall(PetscOptionsHasName(NULL, NULL, "-hard", &flg));
421: if (!flg) {
422: /*
423: Set function evaluation routine and vector.
424: */
425: PetscCall(SNESSetFunction(snes, r, FormFunction1_block, NULL));
427: /*
428: Set Jacobian matrix data structure and Jacobian evaluation routine
429: */
430: PetscCall(SNESSetJacobian(snes, J, J, FormJacobian1_block, NULL));
431: } else {
432: PetscCall(SNESSetFunction(snes, r, FormFunction2_block, NULL));
433: PetscCall(SNESSetJacobian(snes, J, J, FormJacobian2_block, NULL));
434: }
436: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
437: Customize nonlinear solver; set runtime options
438: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
440: /*
441: Set linear solver defaults for this problem. By extracting the
442: KSP, KSP, and PC contexts from the SNES context, we can then
443: directly call any KSP, KSP, and PC routines to set various options.
444: */
445: PetscCall(SNESGetKSP(snes, &ksp));
446: PetscCall(KSPGetPC(ksp, &pc));
447: PetscCall(PCSetType(pc, PCNONE));
448: PetscCall(KSPSetTolerances(ksp, 1.e-4, PETSC_DEFAULT, PETSC_DEFAULT, 20));
450: /*
451: Set SNES/KSP/KSP/PC runtime options, e.g.,
452: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
453: These options will override those specified above as long as
454: SNESSetFromOptions() is called _after_ any other customization
455: routines.
456: */
457: PetscCall(SNESSetFromOptions(snes));
459: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
460: Evaluate initial guess; then solve nonlinear system
461: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
462: if (!flg) {
463: PetscCall(VecSet(x, pfive));
464: } else {
465: Vec *vecs;
466: PetscCall(VecNestGetSubVecs(x, NULL, &vecs));
467: bv = vecs[0];
468: /* PetscCall(VecBlockGetSubVec(x, 0, &bv)); */
469: PetscCall(VecSetValue(bv, 0, 2.0, INSERT_VALUES)); /* xx[0] = 2.0; */
470: PetscCall(VecAssemblyBegin(bv));
471: PetscCall(VecAssemblyEnd(bv));
473: /* PetscCall(VecBlockGetSubVec(x, 1, &bv)); */
474: bv = vecs[1];
475: PetscCall(VecSetValue(bv, 0, 3.0, INSERT_VALUES)); /* xx[1] = 3.0; */
476: PetscCall(VecAssemblyBegin(bv));
477: PetscCall(VecAssemblyEnd(bv));
478: }
479: /*
480: Note: The user should initialize the vector, x, with the initial guess
481: for the nonlinear solver prior to calling SNESSolve(). In particular,
482: to employ an initial guess of zero, the user should explicitly set
483: this vector to zero by calling VecSet().
484: */
485: PetscCall(SNESSolve(snes, NULL, x));
486: PetscCall(SNESGetIterationNumber(snes, &its));
487: if (flg) {
488: Vec f;
489: PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
490: PetscCall(SNESGetFunction(snes, &f, 0, 0));
491: PetscCall(VecView(r, PETSC_VIEWER_STDOUT_WORLD));
492: }
494: PetscCall(PetscPrintf(PETSC_COMM_SELF, "number of SNES iterations = %" PetscInt_FMT "\n\n", its));
496: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
497: Free work space. All PETSc objects should be destroyed when they
498: are no longer needed.
499: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
500: PetscCall(VecDestroy(&x));
501: PetscCall(VecDestroy(&r));
502: PetscCall(MatDestroy(&J));
503: PetscCall(SNESDestroy(&snes));
504: PetscFunctionReturn(PETSC_SUCCESS);
505: }
507: static PetscErrorCode FormFunction1_block(SNES snes, Vec x, Vec f, void *dummy)
508: {
509: Vec *xx, *ff, x1, x2, f1, f2;
510: PetscScalar ff_0, ff_1;
511: PetscScalar xx_0, xx_1;
512: PetscInt index, nb;
514: PetscFunctionBeginUser;
515: /* get blocks for function */
516: PetscCall(VecNestGetSubVecs(f, &nb, &ff));
517: f1 = ff[0];
518: f2 = ff[1];
520: /* get blocks for solution */
521: PetscCall(VecNestGetSubVecs(x, &nb, &xx));
522: x1 = xx[0];
523: x2 = xx[1];
525: /* get solution values */
526: index = 0;
527: PetscCall(VecGetValues(x1, 1, &index, &xx_0));
528: PetscCall(VecGetValues(x2, 1, &index, &xx_1));
530: /* Compute function */
531: ff_0 = xx_0 * xx_0 + xx_0 * xx_1 - 3.0;
532: ff_1 = xx_0 * xx_1 + xx_1 * xx_1 - 6.0;
534: /* set function values */
535: PetscCall(VecSetValue(f1, index, ff_0, INSERT_VALUES));
537: PetscCall(VecSetValue(f2, index, ff_1, INSERT_VALUES));
539: PetscCall(VecAssemblyBegin(f));
540: PetscCall(VecAssemblyEnd(f));
541: PetscFunctionReturn(PETSC_SUCCESS);
542: }
544: static PetscErrorCode FormJacobian1_block(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
545: {
546: Vec *xx, x1, x2;
547: PetscScalar xx_0, xx_1;
548: PetscInt index, nb;
549: PetscScalar A_00, A_01, A_10, A_11;
550: Mat j11, j12, j21, j22;
551: Mat **mats;
553: PetscFunctionBeginUser;
554: /* get blocks for solution */
555: PetscCall(VecNestGetSubVecs(x, &nb, &xx));
556: x1 = xx[0];
557: x2 = xx[1];
559: /* get solution values */
560: index = 0;
561: PetscCall(VecGetValues(x1, 1, &index, &xx_0));
562: PetscCall(VecGetValues(x2, 1, &index, &xx_1));
564: /* get block matrices */
565: PetscCall(MatNestGetSubMats(jac, NULL, NULL, &mats));
566: j11 = mats[0][0];
567: j12 = mats[0][1];
568: j21 = mats[1][0];
569: j22 = mats[1][1];
571: /* compute jacobian entries */
572: A_00 = 2.0 * xx_0 + xx_1;
573: A_01 = xx_0;
574: A_10 = xx_1;
575: A_11 = xx_0 + 2.0 * xx_1;
577: /* set jacobian values */
578: PetscCall(MatSetValue(j11, 0, 0, A_00, INSERT_VALUES));
579: PetscCall(MatSetValue(j12, 0, 0, A_01, INSERT_VALUES));
580: PetscCall(MatSetValue(j21, 0, 0, A_10, INSERT_VALUES));
581: PetscCall(MatSetValue(j22, 0, 0, A_11, INSERT_VALUES));
583: /* Assemble sub matrix */
584: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
585: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
586: PetscFunctionReturn(PETSC_SUCCESS);
587: }
589: static PetscErrorCode FormFunction2_block(SNES snes, Vec x, Vec f, void *dummy)
590: {
591: PetscScalar *ff;
592: const PetscScalar *xx;
594: PetscFunctionBeginUser;
595: /*
596: Get pointers to vector data.
597: - For default PETSc vectors, VecGetArray() returns a pointer to
598: the data array. Otherwise, the routine is implementation dependent.
599: - You MUST call VecRestoreArray() when you no longer need access to
600: the array.
601: */
602: PetscCall(VecGetArrayRead(x, &xx));
603: PetscCall(VecGetArray(f, &ff));
605: /*
606: Compute function
607: */
608: ff[0] = PetscSinScalar(3.0 * xx[0]) + xx[0];
609: ff[1] = xx[1];
611: /*
612: Restore vectors
613: */
614: PetscCall(VecRestoreArrayRead(x, &xx));
615: PetscCall(VecRestoreArray(f, &ff));
616: PetscFunctionReturn(PETSC_SUCCESS);
617: }
619: static PetscErrorCode FormJacobian2_block(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
620: {
621: const PetscScalar *xx;
622: PetscScalar A[4];
623: PetscInt idx[2] = {0, 1};
625: PetscFunctionBeginUser;
626: /*
627: Get pointer to vector data
628: */
629: PetscCall(VecGetArrayRead(x, &xx));
631: /*
632: Compute Jacobian entries and insert into matrix.
633: - Since this is such a small problem, we set all entries for
634: the matrix at once.
635: */
636: A[0] = 3.0 * PetscCosScalar(3.0 * xx[0]) + 1.0;
637: A[1] = 0.0;
638: A[2] = 0.0;
639: A[3] = 1.0;
640: PetscCall(MatSetValues(jac, 2, idx, 2, idx, A, INSERT_VALUES));
642: /*
643: Restore vector
644: */
645: PetscCall(VecRestoreArrayRead(x, &xx));
647: /*
648: Assemble matrix
649: */
650: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
651: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
652: PetscFunctionReturn(PETSC_SUCCESS);
653: }
655: int main(int argc, char **argv)
656: {
657: PetscMPIInt size;
659: PetscFunctionBeginUser;
660: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
661: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
662: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
663: PetscCall(assembled_system());
664: PetscCall(block_system());
665: PetscCall(PetscFinalize());
666: return 0;
667: }
669: /*TEST
671: test:
672: args: -snes_monitor_short
673: requires: !single
675: TEST*/