Actual source code: ex1.c
1: /*$Id: ex1.c,v 1.26 2001/08/07 03:04:16 balay Exp $*/
3: static char help[] = "Newton's method to solve a two-variable system, sequentially.nn";
5: /*T
6: Concepts: SNES^basic uniprocessor example
7: Processors: 1
8: T*/
10: /*
11: Include "petscsnes.h" so that we can use SNES solvers. Note that this
12: file automatically includes:
13: petsc.h - base PETSc routines petscvec.h - vectors
14: petscsys.h - system routines petscmat.h - matrices
15: petscis.h - index sets petscksp.h - Krylov subspace methods
16: petscviewer.h - viewers petscpc.h - preconditioners
17: petscsles.h - linear solvers
18: */
19: #include petscsnes.h
21: /*
22: User-defined routines
23: */
24: extern int FormJacobian1(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
25: extern int FormFunction1(SNES,Vec,Vec,void*);
26: extern int FormJacobian2(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
27: extern int FormFunction2(SNES,Vec,Vec,void*);
29: #undef __FUNCT__
31: int main(int argc,char **argv)
32: {
33: SNES snes; /* nonlinear solver context */
34: SLES sles; /* linear solver context */
35: PC pc; /* preconditioner context */
36: KSP ksp; /* Krylov subspace method context */
37: Vec x,r; /* solution, residual vectors */
38: Mat J; /* Jacobian matrix */
39: int ierr,its,size;
40: PetscScalar pfive = .5,*xx;
41: PetscTruth flg;
43: PetscInitialize(&argc,&argv,(char *)0,help);
44: MPI_Comm_size(PETSC_COMM_WORLD,&size);
45: if (size != 1) SETERRQ(1,"This is a uniprocessor example only!");
47: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
48: Create nonlinear solver context
49: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
51: SNESCreate(PETSC_COMM_WORLD,&snes);
53: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
54: Create matrix and vector data structures; set corresponding routines
55: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
57: /*
58: Create vectors for solution and nonlinear function
59: */
60: VecCreateSeq(PETSC_COMM_SELF,2,&x);
61: VecDuplicate(x,&r);
63: /*
64: Create Jacobian matrix data structure
65: */
66: MatCreate(PETSC_COMM_SELF,PETSC_DECIDE,PETSC_DECIDE,2,2,&J);
67: MatSetFromOptions(J);
69: PetscOptionsHasName(PETSC_NULL,"-hard",&flg);
70: if (!flg) {
71: /*
72: Set function evaluation routine and vector.
73: */
74: SNESSetFunction(snes,r,FormFunction1,PETSC_NULL);
76: /*
77: Set Jacobian matrix data structure and Jacobian evaluation routine
78: */
79: SNESSetJacobian(snes,J,J,FormJacobian1,PETSC_NULL);
80: } else {
81: SNESSetFunction(snes,r,FormFunction2,PETSC_NULL);
82: SNESSetJacobian(snes,J,J,FormJacobian2,PETSC_NULL);
83: }
85: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86: Customize nonlinear solver; set runtime options
87: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
89: /*
90: Set linear solver defaults for this problem. By extracting the
91: SLES, KSP, and PC contexts from the SNES context, we can then
92: directly call any SLES, KSP, and PC routines to set various options.
93: */
94: SNESGetSLES(snes,&sles);
95: SLESGetKSP(sles,&ksp);
96: SLESGetPC(sles,&pc);
97: PCSetType(pc,PCNONE);
98: KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);
100: /*
101: Set SNES/SLES/KSP/PC runtime options, e.g.,
102: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
103: These options will override those specified above as long as
104: SNESSetFromOptions() is called _after_ any other customization
105: routines.
106: */
107: SNESSetFromOptions(snes);
109: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110: Evaluate initial guess; then solve nonlinear system
111: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
112: if (!flg) {
113: VecSet(&pfive,x);
114: } else {
115: VecGetArray(x,&xx);
116: xx[0] = 2.0; xx[1] = 3.0;
117: VecRestoreArray(x,&xx);
118: }
119: /*
120: Note: The user should initialize the vector, x, with the initial guess
121: for the nonlinear solver prior to calling SNESSolve(). In particular,
122: to employ an initial guess of zero, the user should explicitly set
123: this vector to zero by calling VecSet().
124: */
126: SNESSolve(snes,x,&its);
127: if (flg) {
128: Vec f;
129: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
130: SNESGetFunction(snes,&f,0,0);
131: VecView(r,PETSC_VIEWER_STDOUT_WORLD);
132: }
134: PetscPrintf(PETSC_COMM_SELF,"number of Newton iterations = %dnn",its);
136: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137: Free work space. All PETSc objects should be destroyed when they
138: are no longer needed.
139: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
141: VecDestroy(x); VecDestroy(r);
142: MatDestroy(J); SNESDestroy(snes);
144: PetscFinalize();
145: return 0;
146: }
147: /* ------------------------------------------------------------------- */
148: #undef __FUNCT__
150: /*
151: FormFunction1 - Evaluates nonlinear function, F(x).
153: Input Parameters:
154: . snes - the SNES context
155: . x - input vector
156: . dummy - optional user-defined context (not used here)
158: Output Parameter:
159: . f - function vector
160: */
161: int FormFunction1(SNES snes,Vec x,Vec f,void *dummy)
162: {
163: int ierr;
164: PetscScalar *xx,*ff;
166: /*
167: Get pointers to vector data.
168: - For default PETSc vectors, VecGetArray() returns a pointer to
169: the data array. Otherwise, the routine is implementation dependent.
170: - You MUST call VecRestoreArray() when you no longer need access to
171: the array.
172: */
173: VecGetArray(x,&xx);
174: VecGetArray(f,&ff);
176: /*
177: Compute function
178: */
179: ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
180: ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;
182: /*
183: Restore vectors
184: */
185: VecRestoreArray(x,&xx);
186: VecRestoreArray(f,&ff);
188: return 0;
189: }
190: /* ------------------------------------------------------------------- */
191: #undef __FUNCT__
193: /*
194: FormJacobian1 - Evaluates Jacobian matrix.
196: Input Parameters:
197: . snes - the SNES context
198: . x - input vector
199: . dummy - optional user-defined context (not used here)
201: Output Parameters:
202: . jac - Jacobian matrix
203: . B - optionally different preconditioning matrix
204: . flag - flag indicating matrix structure
205: */
206: int FormJacobian1(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
207: {
208: PetscScalar *xx,A[4];
209: int ierr,idx[2] = {0,1};
211: /*
212: Get pointer to vector data
213: */
214: VecGetArray(x,&xx);
216: /*
217: Compute Jacobian entries and insert into matrix.
218: - Since this is such a small problem, we set all entries for
219: the matrix at once.
220: */
221: A[0] = 2.0*xx[0] + xx[1]; A[1] = xx[0];
222: A[2] = xx[1]; A[3] = xx[0] + 2.0*xx[1];
223: MatSetValues(*jac,2,idx,2,idx,A,INSERT_VALUES);
224: *flag = SAME_NONZERO_PATTERN;
226: /*
227: Restore vector
228: */
229: VecRestoreArray(x,&xx);
231: /*
232: Assemble matrix
233: */
234: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
235: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
237: return 0;
238: }
241: /* ------------------------------------------------------------------- */
242: #undef __FUNCT__
244: int FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
245: {
246: int ierr;
247: PetscScalar *xx,*ff;
249: /*
250: Get pointers to vector data.
251: - For default PETSc vectors, VecGetArray() returns a pointer to
252: the data array. Otherwise, the routine is implementation dependent.
253: - You MUST call VecRestoreArray() when you no longer need access to
254: the array.
255: */
256: VecGetArray(x,&xx);
257: VecGetArray(f,&ff);
259: /*
260: Compute function
261: */
262: ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
263: ff[1] = xx[1];
265: /*
266: Restore vectors
267: */
268: VecRestoreArray(x,&xx);
269: VecRestoreArray(f,&ff);
271: return 0;
272: }
273: /* ------------------------------------------------------------------- */
274: #undef __FUNCT__
276: int FormJacobian2(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
277: {
278: PetscScalar *xx,A[4];
279: int ierr,idx[2] = {0,1};
281: /*
282: Get pointer to vector data
283: */
284: VecGetArray(x,&xx);
286: /*
287: Compute Jacobian entries and insert into matrix.
288: - Since this is such a small problem, we set all entries for
289: the matrix at once.
290: */
291: A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
292: A[2] = 0.0; A[3] = 1.0;
293: MatSetValues(*jac,2,idx,2,idx,A,INSERT_VALUES);
294: *flag = SAME_NONZERO_PATTERN;
296: /*
297: Restore vector
298: */
299: VecRestoreArray(x,&xx);
301: /*
302: Assemble matrix
303: */
304: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
305: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
307: return 0;
308: }