Actual source code: ex1.c
2: static char help[] = "Newton's method for a two-variable system, sequential.\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
19: typedef struct {
20: Vec xloc,rloc; /* local solution, residual vectors */
21: VecScatter scatter;
22: } AppCtx;
24: /*
25: User-defined routines
26: */
34: int main(int argc,char **argv)
35: {
36: SNES snes; /* nonlinear solver context */
37: KSP ksp; /* linear solver context */
38: PC pc; /* preconditioner context */
39: Vec x,r; /* solution, residual vectors */
40: Mat J; /* Jacobian matrix */
42: PetscInt its;
43: PetscMPIInt size,rank;
44: PetscScalar pfive = .5,*xx;
45: PetscTruth flg;
46: AppCtx user; /* user-defined work context */
47: IS isglobal,islocal;
49: PetscInitialize(&argc,&argv,(char *)0,help);
50: MPI_Comm_size(PETSC_COMM_WORLD,&size);
51: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
53: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
54: Create nonlinear solver context
55: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
56: SNESCreate(PETSC_COMM_WORLD,&snes);
58: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
59: Create matrix and vector data structures; set corresponding routines
60: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
61: /*
62: Create vectors for solution and nonlinear function
63: */
64: VecCreate(PETSC_COMM_WORLD,&x);
65: VecSetSizes(x,PETSC_DECIDE,2);
66: VecSetFromOptions(x);
67: VecDuplicate(x,&r);
69: if (size > 1){
70: VecCreateSeq(PETSC_COMM_SELF,2,&user.xloc);
71: VecDuplicate(user.xloc,&user.rloc);
73: /* Create the scatter between the global x and local xloc */
74: ISCreateStride(MPI_COMM_SELF,2,0,1,&islocal);
75: ISCreateStride(MPI_COMM_SELF,2,0,1,&isglobal);
76: VecScatterCreate(x,isglobal,user.xloc,islocal,&user.scatter);
77: ISDestroy(isglobal);
78: ISDestroy(islocal);
79: }
81: /*
82: Create Jacobian matrix data structure
83: */
84: MatCreate(PETSC_COMM_WORLD,&J);
85: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
86: MatSetFromOptions(J);
88: PetscOptionsHasName(PETSC_NULL,"-hard",&flg);
89: if (!flg) {
90: /*
91: Set function evaluation routine and vector.
92: */
93: SNESSetFunction(snes,r,FormFunction1,&user);
95: /*
96: Set Jacobian matrix data structure and Jacobian evaluation routine
97: */
98: SNESSetJacobian(snes,J,J,FormJacobian1,PETSC_NULL);
99: } else {
100: if (size != 1) SETERRQ(1,"This case is a uniprocessor example only!");
101: SNESSetFunction(snes,r,FormFunction2,PETSC_NULL);
102: SNESSetJacobian(snes,J,J,FormJacobian2,PETSC_NULL);
103: }
105: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106: Customize nonlinear solver; set runtime options
107: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
108: /*
109: Set linear solver defaults for this problem. By extracting the
110: KSP, KSP, and PC contexts from the SNES context, we can then
111: directly call any KSP, KSP, and PC routines to set various options.
112: */
113: SNESGetKSP(snes,&ksp);
114: KSPGetPC(ksp,&pc);
115: PCSetType(pc,PCNONE);
116: KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);
118: /*
119: Set SNES/KSP/KSP/PC runtime options, e.g.,
120: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
121: These options will override those specified above as long as
122: SNESSetFromOptions() is called _after_ any other customization
123: routines.
124: */
125: SNESSetFromOptions(snes);
127: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: Evaluate initial guess; then solve nonlinear system
129: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130: if (!flg) {
131: VecSet(x,pfive);
132: } else {
133: VecGetArray(x,&xx);
134: xx[0] = 2.0; xx[1] = 3.0;
135: VecRestoreArray(x,&xx);
136: }
137: /*
138: Note: The user should initialize the vector, x, with the initial guess
139: for the nonlinear solver prior to calling SNESSolve(). In particular,
140: to employ an initial guess of zero, the user should explicitly set
141: this vector to zero by calling VecSet().
142: */
144: SNESSolve(snes,PETSC_NULL,x);
145: SNESGetIterationNumber(snes,&its);
146: if (flg) {
147: Vec f;
148: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
149: SNESGetFunction(snes,&f,0,0);
150: VecView(r,PETSC_VIEWER_STDOUT_WORLD);
151: }
153: PetscPrintf(PETSC_COMM_WORLD,"number of Newton iterations = %D\n\n",its);
155: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156: Free work space. All PETSc objects should be destroyed when they
157: are no longer needed.
158: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
160: VecDestroy(x); VecDestroy(r);
161: MatDestroy(J); SNESDestroy(snes);
162: if (size > 1){
163: VecDestroy(user.xloc);
164: VecDestroy(user.rloc);
165: VecScatterDestroy(user.scatter);
166: }
167: PetscFinalize();
168: return 0;
169: }
170: /* ------------------------------------------------------------------- */
173: /*
174: FormFunction1 - Evaluates nonlinear function, F(x).
176: Input Parameters:
177: . snes - the SNES context
178: . x - input vector
179: . ctx - optional user-defined context
181: Output Parameter:
182: . f - function vector
183: */
184: PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ctx)
185: {
187: PetscScalar *xx,*ff;
188: AppCtx *user = (AppCtx*)ctx;
189: Vec xloc=user->xloc,floc=user->rloc;
190: VecScatter scatter=user->scatter;
191: MPI_Comm comm;
192: PetscMPIInt size,rank;
193: PetscInt rstart,rend;
195: PetscObjectGetComm((PetscObject)snes,&comm);
196: MPI_Comm_size(comm,&size);
197: MPI_Comm_rank(comm,&rank);
198: if (size > 1){
199: /*
200: This is a ridiculous case for testing intermidiate steps from sequential
201: code development to parallel implementation.
202: (1) scatter x into a sequetial vector;
203: (2) each process evaluates all values of floc;
204: (3) scatter floc back to the parallel f.
205: */
206: VecScatterBegin(scatter,x,xloc,INSERT_VALUES,SCATTER_FORWARD);
207: VecScatterEnd(scatter,x,xloc,INSERT_VALUES,SCATTER_FORWARD);
209: VecGetOwnershipRange(f,&rstart,&rend);
210: VecGetArray(xloc,&xx);
211: VecGetArray(floc,&ff);
212: ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
213: ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;
214: VecRestoreArray(floc,&ff);
215: VecRestoreArray(xloc,&xx);
217: VecScatterBegin(scatter,floc,f,INSERT_VALUES,SCATTER_REVERSE);
218: VecScatterEnd(scatter,floc,f,INSERT_VALUES,SCATTER_REVERSE);
219: } else {
220: /*
221: Get pointers to vector data.
222: - For default PETSc vectors, VecGetArray() returns a pointer to
223: the data array. Otherwise, the routine is implementation dependent.
224: - You MUST call VecRestoreArray() when you no longer need access to
225: the array.
226: */
227: VecGetArray(x,&xx);
228: VecGetArray(f,&ff);
230: /* Compute function */
231: ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
232: ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;
234: /* Restore vectors */
235: VecRestoreArray(x,&xx);
236: VecRestoreArray(f,&ff);
237: }
238: return 0;
239: }
240: /* ------------------------------------------------------------------- */
243: /*
244: FormJacobian1 - Evaluates Jacobian matrix.
246: Input Parameters:
247: . snes - the SNES context
248: . x - input vector
249: . dummy - optional user-defined context (not used here)
251: Output Parameters:
252: . jac - Jacobian matrix
253: . B - optionally different preconditioning matrix
254: . flag - flag indicating matrix structure
255: */
256: PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
257: {
258: PetscScalar *xx,A[4];
260: PetscInt idx[2] = {0,1};
262: /*
263: Get pointer to vector data
264: */
265: VecGetArray(x,&xx);
267: /*
268: Compute Jacobian entries and insert into matrix.
269: - Since this is such a small problem, we set all entries for
270: the matrix at once.
271: */
272: A[0] = 2.0*xx[0] + xx[1]; A[1] = xx[0];
273: A[2] = xx[1]; A[3] = xx[0] + 2.0*xx[1];
274: MatSetValues(*B,2,idx,2,idx,A,INSERT_VALUES);
275: *flag = SAME_NONZERO_PATTERN;
277: /*
278: Restore vector
279: */
280: VecRestoreArray(x,&xx);
282: /*
283: Assemble matrix
284: */
285: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
286: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
287: if (*jac != *B){
288: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
289: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
290: }
291: return 0;
292: }
294: /* ------------------------------------------------------------------- */
297: PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
298: {
300: PetscScalar *xx,*ff;
302: /*
303: Get pointers to vector data.
304: - For default PETSc vectors, VecGetArray() returns a pointer to
305: the data array. Otherwise, the routine is implementation dependent.
306: - You MUST call VecRestoreArray() when you no longer need access to
307: the array.
308: */
309: VecGetArray(x,&xx);
310: VecGetArray(f,&ff);
312: /*
313: Compute function
314: */
315: ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
316: ff[1] = xx[1];
318: /*
319: Restore vectors
320: */
321: VecRestoreArray(x,&xx);
322: VecRestoreArray(f,&ff);
323: return 0;
324: }
325: /* ------------------------------------------------------------------- */
328: PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
329: {
330: PetscScalar *xx,A[4];
332: PetscInt idx[2] = {0,1};
334: /*
335: Get pointer to vector data
336: */
337: VecGetArray(x,&xx);
339: /*
340: Compute Jacobian entries and insert into matrix.
341: - Since this is such a small problem, we set all entries for
342: the matrix at once.
343: */
344: A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
345: A[2] = 0.0; A[3] = 1.0;
346: MatSetValues(*B,2,idx,2,idx,A,INSERT_VALUES);
347: *flag = SAME_NONZERO_PATTERN;
349: /*
350: Restore vector
351: */
352: VecRestoreArray(x,&xx);
354: /*
355: Assemble matrix
356: */
357: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
358: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
359: if (*jac != *B){
360: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
361: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
362: }
363: return 0;
364: }