Actual source code: ex1.c
2: static char help[] = "Newton's method to solve a two-variable system, sequentially.\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: petsc.h - base PETSc routines petscvec.h - vectors
12: petscsys.h - system routines 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: if (!rank){
154: PetscPrintf(PETSC_COMM_SELF,"number of Newton iterations = %D\n\n",its);
155: }
157: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158: Free work space. All PETSc objects should be destroyed when they
159: are no longer needed.
160: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
162: VecDestroy(x); VecDestroy(r);
163: MatDestroy(J); SNESDestroy(snes);
164: if (size > 1){
165: VecDestroy(user.xloc);
166: VecDestroy(user.rloc);
167: VecScatterDestroy(user.scatter);
168: }
169: PetscFinalize();
170: return 0;
171: }
172: /* ------------------------------------------------------------------- */
175: /*
176: FormFunction1 - Evaluates nonlinear function, F(x).
178: Input Parameters:
179: . snes - the SNES context
180: . x - input vector
181: . ctx - optional user-defined context
183: Output Parameter:
184: . f - function vector
185: */
186: PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ctx)
187: {
189: PetscScalar *xx,*ff;
190: AppCtx *user = (AppCtx*)ctx;
191: Vec xloc=user->xloc,floc=user->rloc;
192: VecScatter scatter=user->scatter;
193: MPI_Comm comm;
194: PetscMPIInt size,rank;
195: PetscInt rstart,rend;
197: PetscObjectGetComm((PetscObject)snes,&comm);
198: MPI_Comm_size(comm,&size);
199: MPI_Comm_rank(comm,&rank);
200: if (size > 1){
201: /*
202: This is a ridiculous case for testing intermidiate steps from sequential code development to parallel implementation.
203: (1) scatter x into a sequetial vector;
204: (2) each process evaluates all values of floc;
205: (3) scatter floc back to the parallel f.
206: */
207: VecScatterBegin(scatter,x,xloc,INSERT_VALUES,SCATTER_FORWARD);
208: VecScatterEnd(scatter,x,xloc,INSERT_VALUES,SCATTER_FORWARD);
210: VecGetOwnershipRange(f,&rstart,&rend);
211: VecGetArray(xloc,&xx);
212: VecGetArray(floc,&ff);
213: ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
214: ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;
215: VecRestoreArray(floc,&xx);
216: VecRestoreArray(xloc,&xx);
218: VecScatterBegin(scatter,floc,f,INSERT_VALUES,SCATTER_REVERSE);
219: VecScatterEnd(scatter,floc,f,INSERT_VALUES,SCATTER_REVERSE);
220: } else {
221: /*
222: Get pointers to vector data.
223: - For default PETSc vectors, VecGetArray() returns a pointer to
224: the data array. Otherwise, the routine is implementation dependent.
225: - You MUST call VecRestoreArray() when you no longer need access to
226: the array.
227: */
228: VecGetArray(x,&xx);
229: VecGetArray(f,&ff);
231: /* Compute function */
232: ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
233: ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;
235: /* Restore vectors */
236: VecRestoreArray(x,&xx);
237: VecRestoreArray(f,&ff);
238: }
239: return 0;
240: }
241: /* ------------------------------------------------------------------- */
244: /*
245: FormJacobian1 - Evaluates Jacobian matrix.
247: Input Parameters:
248: . snes - the SNES context
249: . x - input vector
250: . dummy - optional user-defined context (not used here)
252: Output Parameters:
253: . jac - Jacobian matrix
254: . B - optionally different preconditioning matrix
255: . flag - flag indicating matrix structure
256: */
257: PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
258: {
259: PetscScalar *xx,A[4];
261: PetscInt idx[2] = {0,1};
263: /*
264: Get pointer to vector data
265: */
266: VecGetArray(x,&xx);
268: /*
269: Compute Jacobian entries and insert into matrix.
270: - Since this is such a small problem, we set all entries for
271: the matrix at once.
272: */
273: A[0] = 2.0*xx[0] + xx[1]; A[1] = xx[0];
274: A[2] = xx[1]; A[3] = xx[0] + 2.0*xx[1];
275: MatSetValues(*B,2,idx,2,idx,A,INSERT_VALUES);
276: *flag = SAME_NONZERO_PATTERN;
278: /*
279: Restore vector
280: */
281: VecRestoreArray(x,&xx);
283: /*
284: Assemble matrix
285: */
286: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
287: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
288: if (*jac != *B){
289: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
290: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
291: }
292: return 0;
293: }
295: /* ------------------------------------------------------------------- */
298: PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
299: {
301: PetscScalar *xx,*ff;
303: /*
304: Get pointers to vector data.
305: - For default PETSc vectors, VecGetArray() returns a pointer to
306: the data array. Otherwise, the routine is implementation dependent.
307: - You MUST call VecRestoreArray() when you no longer need access to
308: the array.
309: */
310: VecGetArray(x,&xx);
311: VecGetArray(f,&ff);
313: /*
314: Compute function
315: */
316: ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
317: ff[1] = xx[1];
319: /*
320: Restore vectors
321: */
322: VecRestoreArray(x,&xx);
323: VecRestoreArray(f,&ff);
324: return 0;
325: }
326: /* ------------------------------------------------------------------- */
329: PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
330: {
331: PetscScalar *xx,A[4];
333: PetscInt idx[2] = {0,1};
335: /*
336: Get pointer to vector data
337: */
338: VecGetArray(x,&xx);
340: /*
341: Compute Jacobian entries and insert into matrix.
342: - Since this is such a small problem, we set all entries for
343: the matrix at once.
344: */
345: A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
346: A[2] = 0.0; A[3] = 1.0;
347: MatSetValues(*B,2,idx,2,idx,A,INSERT_VALUES);
348: *flag = SAME_NONZERO_PATTERN;
350: /*
351: Restore vector
352: */
353: VecRestoreArray(x,&xx);
355: /*
356: Assemble matrix
357: */
358: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
359: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
360: if (*jac != *B){
361: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
362: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
363: }
364: return 0;
365: }