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;
33: PetscInitialize(&argc,&argv,(char *)0,help);
34:
35: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
36: Create nonlinear solver context
37: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
38: SNESCreate(PETSC_COMM_WORLD,&snes);
40: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
41: Create matrix and vector data structures; set corresponding routines
42: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
43: /*
44: Create vectors for solution and nonlinear function
45: */
46: VecCreate(PETSC_COMM_WORLD,&x);
47: VecSetSizes(x,PETSC_DECIDE,2);
48: VecSetFromOptions(x);
49: VecDuplicate(x,&r);
51: /*
52: Create Jacobian matrix data structure
53: */
54: MatCreate(PETSC_COMM_WORLD,&J);
55: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
56: MatSetFromOptions(J);
58: /*
59: Set function evaluation routine and vector.
60: */
61: SNESSetFunction(snes,r,FormFunction1,PETSC_NULL);
63: /*
64: Set Jacobian matrix data structure and Jacobian evaluation routine
65: */
66: SNESSetJacobian(snes,J,J,FormJacobian1,PETSC_NULL);
68: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69: Customize nonlinear solver; set runtime options
70: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
71: SNESSetFromOptions(snes);
73: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
74: Evaluate initial guess; then solve nonlinear system
75: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
76: VecGetArray(x,&xx);
77: xx[0] = -1.2; xx[1] = 1.0;
78: VecRestoreArray(x,&xx);
80: /*
81: Note: The user should initialize the vector, x, with the initial guess
82: for the nonlinear solver prior to calling SNESSolve(). In particular,
83: to employ an initial guess of zero, the user should explicitly set
84: this vector to zero by calling VecSet().
85: */
87: SNESSolve(snes,PETSC_NULL,x);
88: VecView(x,0);
89: SNESGetIterationNumber(snes,&its);
90: PetscPrintf(PETSC_COMM_SELF,"number of Newton iterations = %D\n\n",its);
92: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
93: Free work space. All PETSc objects should be destroyed when they
94: are no longer needed.
95: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
97: VecDestroy(x); VecDestroy(r);
98: MatDestroy(J); SNESDestroy(snes);
99: PetscFinalize();
100: return 0;
101: }
102: /* ------------------------------------------------------------------- */
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: {
119: PetscScalar *xx,*ff;
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: VecGetArray(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: VecRestoreArray(x,&xx);
137: VecRestoreArray(f,&ff);
138: return 0;
139: }
140: /* ------------------------------------------------------------------- */
143: /*
144: FormJacobian1 - Evaluates Jacobian matrix.
146: Input Parameters:
147: . snes - the SNES context
148: . x - input vector
149: . dummy - optional user-defined context (not used here)
151: Output Parameters:
152: . jac - Jacobian matrix
153: . B - optionally different preconditioning matrix
154: . flag - flag indicating matrix structure
155: */
156: PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
157: {
158: PetscScalar *xx,A[4];
160: PetscInt idx[2] = {0,1};
162: /*
163: Get pointer to vector data
164: */
165: VecGetArray(x,&xx);
167: /*
168: Compute Jacobian entries and insert into matrix.
169: - Since this is such a small problem, we set all entries for
170: the matrix at once.
171: */
172: A[0] = 2.0 + 1200.0*xx[0]*xx[0] - 400.0*xx[1];
173: A[1] = -400.0*xx[0];
174: A[2] = -400.0*xx[0];
175: A[3] = 200;
176: MatSetValues(*B,2,idx,2,idx,A,INSERT_VALUES);
177: *flag = SAME_NONZERO_PATTERN;
179: /*
180: Restore vector
181: */
182: VecRestoreArray(x,&xx);
184: /*
185: Assemble matrix
186: */
187: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
188: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
189: if (*jac != *B){
190: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
191: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
192: }
193: return 0;
194: }