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: /*T
5: Concepts: SNES^basic example
6: T*/
8: /*
9: Include "petscsnes.h" so that we can use SNES solvers. Note that this
10: file automatically includes:
11: petscsys.h - base PETSc routines petscvec.h - vectors
12: petscmat.h - matrices
13: petscis.h - index sets petscksp.h - Krylov subspace methods
14: petscviewer.h - viewers petscpc.h - preconditioners
15: petscksp.h - linear solvers
16: */
17: #include <petscsnes.h>
24: int main(int argc,char **argv)
25: {
26: SNES snes; /* nonlinear solver context */
27: Vec x,r; /* solution, residual vectors */
28: Mat J; /* Jacobian matrix */
30: PetscInt its;
31: PetscScalar *xx;
32: SNESConvergedReason reason;
34: PetscInitialize(&argc,&argv,(char *)0,help);
35:
36: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
37: Create nonlinear solver context
38: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
39: SNESCreate(PETSC_COMM_WORLD,&snes);
41: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
42: Create matrix and vector data structures; set corresponding routines
43: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
44: /*
45: Create vectors for solution and nonlinear function
46: */
47: VecCreate(PETSC_COMM_WORLD,&x);
48: VecSetSizes(x,PETSC_DECIDE,2);
49: VecSetFromOptions(x);
50: VecDuplicate(x,&r);
52: /*
53: Create Jacobian matrix data structure
54: */
55: MatCreate(PETSC_COMM_WORLD,&J);
56: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
57: MatSetFromOptions(J);
59: /*
60: Set function evaluation routine and vector.
61: */
62: SNESSetFunction(snes,r,FormFunction1,PETSC_NULL);
64: /*
65: Set Jacobian matrix data structure and Jacobian evaluation routine
66: */
67: SNESSetJacobian(snes,J,J,FormJacobian1,PETSC_NULL);
69: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
70: Customize nonlinear solver; set runtime options
71: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
72: SNESSetFromOptions(snes);
74: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
75: Evaluate initial guess; then solve nonlinear system
76: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
77: VecGetArray(x,&xx);
78: xx[0] = -1.2; xx[1] = 1.0;
79: VecRestoreArray(x,&xx);
81: /*
82: Note: The user should initialize the vector, x, with the initial guess
83: for the nonlinear solver prior to calling SNESSolve(). In particular,
84: to employ an initial guess of zero, the user should explicitly set
85: this vector to zero by calling VecSet().
86: */
88: SNESSolve(snes,PETSC_NULL,x);
89: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
90: SNESGetIterationNumber(snes,&its);
91: SNESGetConvergedReason(snes,&reason);
92: /*
93: Some systems computes a residual that is identically zero, thus converging
94: due to FNORM_ABS, others converge due to FNORM_RELATIVE. Here, we only
95: report the reason if the iteration did not converge so that the tests are
96: reproducible.
97: */
98: PetscPrintf(PETSC_COMM_WORLD,"%s number of Newton iterations = %D\n\n",reason>0?"CONVERGED":SNESConvergedReasons[reason],its);
100: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: Free work space. All PETSc objects should be destroyed when they
102: are no longer needed.
103: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
105: VecDestroy(&x); VecDestroy(&r);
106: MatDestroy(&J); SNESDestroy(&snes);
107: PetscFinalize();
108: return 0;
109: }
110: /* ------------------------------------------------------------------- */
113: /*
114: FormFunction1 - Evaluates nonlinear function, F(x).
116: Input Parameters:
117: . snes - the SNES context
118: . x - input vector
119: . ctx - optional user-defined context
121: Output Parameter:
122: . f - function vector
123: */
124: PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ctx)
125: {
127: PetscScalar *xx,*ff;
129: /*
130: Get pointers to vector data.
131: - For default PETSc vectors, VecGetArray() returns a pointer to
132: the data array. Otherwise, the routine is implementation dependent.
133: - You MUST call VecRestoreArray() when you no longer need access to
134: the array.
135: */
136: VecGetArray(x,&xx);
137: VecGetArray(f,&ff);
139: /* Compute function */
140: ff[0] = -2.0 + 2.0*xx[0] + 400.0*xx[0]*xx[0]*xx[0] - 400.0*xx[0]*xx[1];
141: ff[1] = -200.0*xx[0]*xx[0] + 200.0*xx[1];
143: /* Restore vectors */
144: VecRestoreArray(x,&xx);
145: VecRestoreArray(f,&ff);
146: return 0;
147: }
148: /* ------------------------------------------------------------------- */
151: /*
152: FormJacobian1 - Evaluates Jacobian matrix.
154: Input Parameters:
155: . snes - the SNES context
156: . x - input vector
157: . dummy - optional user-defined context (not used here)
159: Output Parameters:
160: . jac - Jacobian matrix
161: . B - optionally different preconditioning matrix
162: . flag - flag indicating matrix structure
163: */
164: PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
165: {
166: PetscScalar *xx,A[4];
168: PetscInt idx[2] = {0,1};
170: /*
171: Get pointer to vector data
172: */
173: VecGetArray(x,&xx);
175: /*
176: Compute Jacobian entries and insert into matrix.
177: - Since this is such a small problem, we set all entries for
178: the matrix at once.
179: */
180: A[0] = 2.0 + 1200.0*xx[0]*xx[0] - 400.0*xx[1];
181: A[1] = -400.0*xx[0];
182: A[2] = -400.0*xx[0];
183: A[3] = 200;
184: MatSetValues(*B,2,idx,2,idx,A,INSERT_VALUES);
185: *flag = SAME_NONZERO_PATTERN;
187: /*
188: Restore vector
189: */
190: VecRestoreArray(x,&xx);
192: /*
193: Assemble matrix
194: */
195: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
196: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
197: if (*jac != *B){
198: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
199: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
200: }
201: return 0;
202: }