Actual source code: ex42.c
2: static char help[] = "Newton's method to solve a two-variable system that comes from the Rosenbrock function.\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: 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: extern PetscErrorCode FormJacobian1(SNES,Vec,Mat,Mat,void*);
16: extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);
18: int main(int argc,char **argv)
19: {
20: SNES snes; /* nonlinear solver context */
21: Vec x,r; /* solution, residual vectors */
22: Mat J; /* Jacobian matrix */
23: PetscInt its;
24: PetscScalar *xx;
25: SNESConvergedReason reason;
27: PetscInitialize(&argc,&argv,(char*)0,help);
29: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
30: Create nonlinear solver context
31: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
32: SNESCreate(PETSC_COMM_WORLD,&snes);
34: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
35: Create matrix and vector data structures; set corresponding routines
36: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
37: /*
38: Create vectors for solution and nonlinear function
39: */
40: VecCreate(PETSC_COMM_WORLD,&x);
41: VecSetSizes(x,PETSC_DECIDE,2);
42: VecSetFromOptions(x);
43: VecDuplicate(x,&r);
45: /*
46: Create Jacobian matrix data structure
47: */
48: MatCreate(PETSC_COMM_WORLD,&J);
49: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
50: MatSetFromOptions(J);
51: MatSetUp(J);
53: /*
54: Set function evaluation routine and vector.
55: */
56: SNESSetFunction(snes,r,FormFunction1,NULL);
58: /*
59: Set Jacobian matrix data structure and Jacobian evaluation routine
60: */
61: SNESSetJacobian(snes,J,J,FormJacobian1,NULL);
63: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
64: Customize nonlinear solver; set runtime options
65: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
66: SNESSetFromOptions(snes);
68: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69: Evaluate initial guess; then solve nonlinear system
70: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
71: VecGetArray(x,&xx);
72: xx[0] = -1.2; xx[1] = 1.0;
73: VecRestoreArray(x,&xx);
75: /*
76: Note: The user should initialize the vector, x, with the initial guess
77: for the nonlinear solver prior to calling SNESSolve(). In particular,
78: to employ an initial guess of zero, the user should explicitly set
79: this vector to zero by calling VecSet().
80: */
82: SNESSolve(snes,NULL,x);
83: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
84: SNESGetIterationNumber(snes,&its);
85: SNESGetConvergedReason(snes,&reason);
86: /*
87: Some systems computes a residual that is identically zero, thus converging
88: due to FNORM_ABS, others converge due to FNORM_RELATIVE. Here, we only
89: report the reason if the iteration did not converge so that the tests are
90: reproducible.
91: */
92: PetscPrintf(PETSC_COMM_WORLD,"%s number of SNES iterations = %D\n",reason>0 ? "CONVERGED" : SNESConvergedReasons[reason],its);
94: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: Free work space. All PETSc objects should be destroyed when they
96: are no longer needed.
97: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99: VecDestroy(&x)); PetscCall(VecDestroy(&r);
100: MatDestroy(&J)); PetscCall(SNESDestroy(&snes);
101: PetscFinalize();
102: return 0;
103: }
104: /* ------------------------------------------------------------------- */
105: /*
106: FormFunction1 - Evaluates nonlinear function, F(x).
108: Input Parameters:
109: . snes - the SNES context
110: . x - input vector
111: . ctx - optional user-defined context
113: Output Parameter:
114: . f - function vector
115: */
116: PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ctx)
117: {
118: PetscScalar *ff;
119: const PetscScalar *xx;
121: /*
122: Get pointers to vector data.
123: - For default PETSc vectors, VecGetArray() returns a pointer to
124: the data array. Otherwise, the routine is implementation dependent.
125: - You MUST call VecRestoreArray() when you no longer need access to
126: the array.
127: */
128: VecGetArrayRead(x,&xx);
129: VecGetArray(f,&ff);
131: /* Compute function */
132: ff[0] = -2.0 + 2.0*xx[0] + 400.0*xx[0]*xx[0]*xx[0] - 400.0*xx[0]*xx[1];
133: ff[1] = -200.0*xx[0]*xx[0] + 200.0*xx[1];
135: /* Restore vectors */
136: VecRestoreArrayRead(x,&xx);
137: VecRestoreArray(f,&ff);
138: return 0;
139: }
140: /* ------------------------------------------------------------------- */
141: /*
142: FormJacobian1 - Evaluates Jacobian matrix.
144: Input Parameters:
145: . snes - the SNES context
146: . x - input vector
147: . dummy - optional user-defined context (not used here)
149: Output Parameters:
150: . jac - Jacobian matrix
151: . B - optionally different preconditioning matrix
152: . flag - flag indicating matrix structure
153: */
154: PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
155: {
156: const PetscScalar *xx;
157: PetscScalar A[4];
158: PetscInt idx[2] = {0,1};
160: /*
161: Get pointer to vector data
162: */
163: VecGetArrayRead(x,&xx);
165: /*
166: Compute Jacobian entries and insert into matrix.
167: - Since this is such a small problem, we set all entries for
168: the matrix at once.
169: */
170: A[0] = 2.0 + 1200.0*xx[0]*xx[0] - 400.0*xx[1];
171: A[1] = -400.0*xx[0];
172: A[2] = -400.0*xx[0];
173: A[3] = 200;
174: MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES);
176: /*
177: Restore vector
178: */
179: VecRestoreArrayRead(x,&xx);
181: /*
182: Assemble matrix
183: */
184: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
185: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
186: if (jac != B) {
187: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
188: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
189: }
190: return 0;
191: }
193: /*TEST
195: test:
196: suffix: 1
197: args: -snes_monitor_short -snes_max_it 1000
198: requires: !single
200: test:
201: suffix: 2
202: args: -snes_monitor_short -snes_max_it 1000 -snes_type newtontrdc -snes_trdc_use_cauchy false
203: requires: !single
205: TEST*/