Actual source code: ex43.c
petsc-3.3-p7 2013-05-11
2: static char help[] = "Newton's method to solve a many-variable system that comes from the 2 variable Rosenbrock function + trivial.\n\n";
4: /*
6: ./ex43 -snes_monitor_range -snes_max_it 1000 -snes_rtol 1.e-14 -n 10 -snes_converged_reason -sub_snes_monito -sub_snes_mf -sub_snes_converged_reason -sub_snes_rtol 1.e-10 -sub_snes_max_it 1000 -sub_snes_monitor
8: Accelerates Newton's method by solving a small problem defined by those elements with large residual plus one level of overlap
10: This is a toy code for playing around
12: Counts residual entries as small if they are less then .2 times the maximum
13: Decides to solve a reduced problem if the number of large entries is less than 20 percent of all entries (and this has been true for criteria_reduce iterations)
14: */
15: #include "ex43-44.h"
18: extern PetscErrorCode FormJacobian1(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
19: extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);
21: typedef struct {
22: PetscInt n,p;
23: } Ctx;
27: int main(int argc,char **argv)
28: {
29: SNES snes; /* nonlinear solver context */
30: Vec x,r; /* solution, residual vectors */
31: Mat J; /* Jacobian matrix */
32: PetscErrorCode ierr;
33: PetscScalar *xx;
34: PetscInt i,max_snes_solves = 20,snes_steps_per_solve = 2 ,criteria_reduce = 1;
35: Ctx ctx;
36: SNESConvergedReason reason;
38: PetscInitialize(&argc,&argv,(char *)0,help);
39: ctx.n = 0;
40: PetscOptionsGetInt(PETSC_NULL,"-n",&ctx.n,PETSC_NULL);
41: ctx.p = 0;
42: PetscOptionsGetInt(PETSC_NULL,"-p",&ctx.p,PETSC_NULL);
43: PetscOptionsGetInt(PETSC_NULL,"-max_snes_solves",&max_snes_solves,PETSC_NULL);
44: PetscOptionsGetInt(PETSC_NULL,"-snes_steps_per_solve",&snes_steps_per_solve,PETSC_NULL);
45: PetscOptionsGetInt(PETSC_NULL,"-criteria_reduce",&criteria_reduce,PETSC_NULL);
47: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
48: Create nonlinear solver context
49: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
50: SNESCreate(PETSC_COMM_WORLD,&snes);
52: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
53: Create matrix and vector data structures; set corresponding routines
54: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
55: /*
56: Create vectors for solution and nonlinear function
57: */
58: VecCreate(PETSC_COMM_WORLD,&x);
59: VecSetSizes(x,PETSC_DECIDE,2+ctx.n+ctx.p);
60: VecSetFromOptions(x);
61: VecDuplicate(x,&r);
63: /*
64: Create Jacobian matrix data structure
65: */
66: MatCreate(PETSC_COMM_WORLD,&J);
67: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2+ctx.p+ctx.n,2+ctx.p+ctx.n);
68: MatSetFromOptions(J);
70: /*
71: Set function evaluation routine and vector.
72: */
73: SNESSetFunction(snes,r,FormFunction1,(void*)&ctx);
75: /*
76: Set Jacobian matrix data structure and Jacobian evaluation routine
77: */
78: SNESSetJacobian(snes,J,J,FormJacobian1,(void*)&ctx);
80: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
81: Customize nonlinear solver; set runtime options
82: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
83: SNESSetFromOptions(snes);
85: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86: Evaluate initial guess; then solve nonlinear system
87: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
88: VecSet(x,0.0);
89: VecGetArray(x,&xx);
90: xx[0] = -1.2;
91: for (i=1; i<ctx.p+2; i++) {
92: xx[i] = 1.0;
93: }
94: VecRestoreArray(x,&xx);
96: /*
97: Note: The user should initialize the vector, x, with the initial guess
98: for the nonlinear solver prior to calling SNESSolve(). In particular,
99: to employ an initial guess of zero, the user should explicitly set
100: this vector to zero by calling VecSet().
101: */
103: SNESMonitorSet(snes,MonitorRange,0,0);
104: SNESSetTolerances(snes,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,snes_steps_per_solve,PETSC_DEFAULT);
105: for (i=0; i<max_snes_solves; i++) {
106: SNESSolve(snes,PETSC_NULL,x);
107: SNESGetConvergedReason(snes,&reason);
108: if (reason && reason != SNES_DIVERGED_MAX_IT) break;
109: if (CountGood > criteria_reduce) {
110: SolveSubproblem(snes);
111: CountGood = 0;
112: }
113: }
115: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: Free work space. All PETSc objects should be destroyed when they
117: are no longer needed.
118: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120: VecDestroy(&x); VecDestroy(&r);
121: MatDestroy(&J); SNESDestroy(&snes);
122: PetscFinalize();
123: return 0;
124: }
125: /* ------------------------------------------------------------------- */
128: /*
129: FormFunction1 - Evaluates nonlinear function, F(x).
131: Input Parameters:
132: . snes - the SNES context
133: . x - input vector
134: . ctx - optional user-defined context
136: Output Parameter:
137: . f - function vector
138: */
139: PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ictx)
140: {
142: PetscScalar *xx,*ff;
143: PetscInt i;
144: Ctx *ctx = (Ctx*)ictx;
146: /*
147: Get pointers to vector data.
148: - For default PETSc vectors, VecGetArray() returns a pointer to
149: the data array. Otherwise, the routine is implementation dependent.
150: - You MUST call VecRestoreArray() when you no longer need access to
151: the array.
152: */
153: VecGetArray(x,&xx);
154: VecGetArray(f,&ff);
156: /* Compute function */
157: ff[0] = -2.0 + 2.0*xx[0] + 400.0*xx[0]*xx[0]*xx[0] - 400.0*xx[0]*xx[1];
158: for (i=1; i<1+ctx->p; i++) {
159: ff[i] = -2.0 + 2.0*xx[i] + 400.0*xx[i]*xx[i]*xx[i] - 400.0*xx[i]*xx[i+1] + 200.0*(xx[i] - xx[i-1]*xx[i-1]);
160: }
161: ff[ctx->p+1] = -200.0*xx[ctx->p]*xx[ctx->p] + 200.0*xx[ctx->p+1];
162: for (i=ctx->p+2; i<2+ctx->p+ctx->n; i++) {
163: ff[i] = xx[i] - xx[0] + .7*xx[1] - .2*xx[i-1]*xx[i-1];
164: }
166: /* Restore vectors */
167: VecRestoreArray(x,&xx);
168: VecRestoreArray(f,&ff);
169: return 0;
170: }
171: /* ------------------------------------------------------------------- */
174: /*
175: FormJacobian1 - Evaluates Jacobian matrix.
177: Input Parameters:
178: . snes - the SNES context
179: . x - input vector
180: . dummy - optional user-defined context (not used here)
182: Output Parameters:
183: . jac - Jacobian matrix
184: . B - optionally different preconditioning matrix
185: . flag - flag indicating matrix structure
186: */
187: PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *ictx)
188: {
189: PetscScalar *xx;
191: PetscInt i;
192: Ctx *ctx = (Ctx*)ictx;
194: MatZeroEntries(*B);
195: /*
196: Get pointer to vector data
197: */
198: VecGetArray(x,&xx);
200: /*
201: Compute Jacobian entries and insert into matrix.
202: - Since this is such a small problem, we set all entries for
203: the matrix at once.
204: */
205: MatSetValue(*B,0,0, 2.0 + 1200.0*xx[0]*xx[0] - 400.0*xx[1],ADD_VALUES);
206: MatSetValue(*B,0,1,-400.0*xx[0] ,ADD_VALUES);
208: for (i=1; i<ctx->p+1; i++) {
209: MatSetValue(*B,i,i-1, -400.0*xx[i-1] ,ADD_VALUES);
210: MatSetValue(*B,i,i, 2.0 + 1200.0*xx[i]*xx[i] - 400.0*xx[i+1] + 200.0,ADD_VALUES);
211: MatSetValue(*B,i,i+1,-400.0*xx[i] ,ADD_VALUES);
212: }
214: MatSetValue(*B,ctx->p+1,ctx->p, -400.0*xx[ctx->p] ,ADD_VALUES);
215: MatSetValue(*B,ctx->p+1,ctx->p+1,200 ,ADD_VALUES);
217: *flag = SAME_NONZERO_PATTERN;
219: for (i=ctx->p+2; i<2+ctx->p+ctx->n; i++) {
220: MatSetValue(*B,i,i,1.0,ADD_VALUES);
221: MatSetValue(*B,i,0,-1.0,ADD_VALUES);
222: MatSetValue(*B,i,1,.7,ADD_VALUES);
223: MatSetValue(*B,i,i-1,-.4*xx[i-1],ADD_VALUES);
224: }
225: /*
226: Restore vector
227: */
228: VecRestoreArray(x,&xx);
230: /*
231: Assemble matrix
232: */
233: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
234: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
235: if (*jac != *B){
236: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
237: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
238: }
239: return 0;
240: }