Actual source code: rosenbrock3.c
1: /* Program usage: mpiexec -n 1 rosenbrock2 [-help] [all TAO options] */
3: /* Include "petsctao.h" so we can use TAO solvers. */
4: #include <petsctao.h>
6: static char help[] = "This example demonstrates use of the TAO package to \n\
7: solve an unconstrained minimization problem on a single processor. We \n\
8: minimize the extended Rosenbrock function: \n\
9: sum_{i=0}^{n/2-1} (alpha*(x_{2i+1}-x_{2i}^2)^2 + (1-x_{2i})^2) \n\
10: or the chained Rosenbrock function:\n\
11: sum_{i=0}^{n-1} alpha*(x_{i+1} - x_i^2)^2 + (1 - x_i)^2\n";
13: /*
14: User-defined application context - contains data needed by the
15: application-provided call-back routines that evaluate the function,
16: gradient, and hessian.
17: */
18: typedef struct {
19: PetscInt n; /* dimension */
20: PetscReal alpha; /* condition parameter */
21: PetscBool chained;
22: } AppCtx;
24: /* -------------- User-defined routines ---------- */
25: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*);
26: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*);
28: int main(int argc,char **argv)
29: {
30: PetscReal zero=0.0;
31: Vec x; /* solution vector */
32: Mat H;
33: Tao tao; /* Tao solver context */
34: PetscBool flg, test_lmvm = PETSC_FALSE;
35: PetscMPIInt size; /* number of processes running */
36: AppCtx user; /* user-defined application context */
37: TaoConvergedReason reason;
38: PetscInt its, recycled_its=0, oneshot_its=0;
40: /* Initialize TAO and PETSc */
41: PetscInitialize(&argc,&argv,(char*)0,help);
42: MPI_Comm_size(PETSC_COMM_WORLD,&size);
45: /* Initialize problem parameters */
46: user.n = 2; user.alpha = 99.0; user.chained = PETSC_FALSE;
47: /* Check for command line arguments to override defaults */
48: PetscOptionsGetInt(NULL,NULL,"-n",&user.n,&flg);
49: PetscOptionsGetReal(NULL,NULL,"-alpha",&user.alpha,&flg);
50: PetscOptionsGetBool(NULL,NULL,"-chained",&user.chained,&flg);
51: PetscOptionsGetBool(NULL,NULL,"-test_lmvm",&test_lmvm,&flg);
53: /* Allocate vectors for the solution and gradient */
54: VecCreateSeq(PETSC_COMM_SELF,user.n,&x);
55: MatCreateSeqBAIJ(PETSC_COMM_SELF,2,user.n,user.n,1,NULL,&H);
57: /* The TAO code begins here */
59: /* Create TAO solver with desired solution method */
60: TaoCreate(PETSC_COMM_SELF,&tao);
61: TaoSetType(tao,TAOBQNLS);
63: /* Set solution vec and an initial guess */
64: VecSet(x, zero);
65: TaoSetSolution(tao,x);
67: /* Set routines for function, gradient, hessian evaluation */
68: TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,&user);
69: TaoSetHessian(tao,H,H,FormHessian,&user);
71: /* Check for TAO command line options */
72: TaoSetFromOptions(tao);
74: /* Solve the problem */
75: TaoSetTolerances(tao, 1.e-5, 0.0, 0.0);
76: TaoSetMaximumIterations(tao, 5);
77: TaoSetRecycleHistory(tao, PETSC_TRUE);
78: reason = TAO_CONTINUE_ITERATING;
79: flg = PETSC_FALSE;
80: TaoGetRecycleHistory(tao, &flg);
81: if (flg) PetscPrintf(PETSC_COMM_SELF, "Recycle: enabled\n");
82: while (reason != TAO_CONVERGED_GATOL) {
83: TaoSolve(tao);
84: TaoGetConvergedReason(tao, &reason);
85: TaoGetIterationNumber(tao, &its);
86: recycled_its += its;
87: PetscPrintf(PETSC_COMM_SELF, "-----------------------\n");
88: }
90: /* Disable recycling and solve again! */
91: TaoSetMaximumIterations(tao, 100);
92: TaoSetRecycleHistory(tao, PETSC_FALSE);
93: VecSet(x, zero);
94: TaoGetRecycleHistory(tao, &flg);
95: if (!flg) PetscPrintf(PETSC_COMM_SELF, "Recycle: disabled\n");
96: TaoSolve(tao);
97: TaoGetConvergedReason(tao, &reason);
99: TaoGetIterationNumber(tao, &oneshot_its);
100: PetscPrintf(PETSC_COMM_SELF, "-----------------------\n");
101: PetscPrintf(PETSC_COMM_SELF, "recycled its: %D | oneshot its: %D\n", recycled_its, oneshot_its);
104: TaoDestroy(&tao);
105: VecDestroy(&x);
106: MatDestroy(&H);
108: PetscFinalize();
109: return 0;
110: }
112: /* -------------------------------------------------------------------- */
113: /*
114: FormFunctionGradient - Evaluates the function, f(X), and gradient, G(X).
116: Input Parameters:
117: . tao - the Tao context
118: . X - input vector
119: . ptr - optional user-defined context, as set by TaoSetFunctionGradient()
121: Output Parameters:
122: . G - vector containing the newly evaluated gradient
123: . f - function value
125: Note:
126: Some optimization methods ask for the function and the gradient evaluation
127: at the same time. Evaluating both at once may be more efficient than
128: evaluating each separately.
129: */
130: PetscErrorCode FormFunctionGradient(Tao tao,Vec X,PetscReal *f, Vec G,void *ptr)
131: {
132: AppCtx *user = (AppCtx *) ptr;
133: PetscInt i,nn=user->n/2;
134: PetscReal ff=0,t1,t2,alpha=user->alpha;
135: PetscScalar *g;
136: const PetscScalar *x;
139: /* Get pointers to vector data */
140: VecGetArrayRead(X,&x);
141: VecGetArrayWrite(G,&g);
143: /* Compute G(X) */
144: if (user->chained) {
145: g[0] = 0;
146: for (i=0; i<user->n-1; i++) {
147: t1 = x[i+1] - x[i]*x[i];
148: ff += PetscSqr(1 - x[i]) + alpha*t1*t1;
149: g[i] += -2*(1 - x[i]) + 2*alpha*t1*(-2*x[i]);
150: g[i+1] = 2*alpha*t1;
151: }
152: } else {
153: for (i=0; i<nn; i++) {
154: t1 = x[2*i+1]-x[2*i]*x[2*i]; t2= 1-x[2*i];
155: ff += alpha*t1*t1 + t2*t2;
156: g[2*i] = -4*alpha*t1*x[2*i]-2.0*t2;
157: g[2*i+1] = 2*alpha*t1;
158: }
159: }
161: /* Restore vectors */
162: VecRestoreArrayRead(X,&x);
163: VecRestoreArrayWrite(G,&g);
164: *f = ff;
166: PetscLogFlops(15.0*nn);
167: return 0;
168: }
170: /* ------------------------------------------------------------------- */
171: /*
172: FormHessian - Evaluates Hessian matrix.
174: Input Parameters:
175: . tao - the Tao context
176: . x - input vector
177: . ptr - optional user-defined context, as set by TaoSetHessian()
179: Output Parameters:
180: . H - Hessian matrix
182: Note: Providing the Hessian may not be necessary. Only some solvers
183: require this matrix.
184: */
185: PetscErrorCode FormHessian(Tao tao,Vec X,Mat H, Mat Hpre, void *ptr)
186: {
187: AppCtx *user = (AppCtx*)ptr;
188: PetscInt i, ind[2];
189: PetscReal alpha=user->alpha;
190: PetscReal v[2][2];
191: const PetscScalar *x;
192: PetscBool assembled;
195: /* Zero existing matrix entries */
196: MatAssembled(H,&assembled);
197: if (assembled || user->chained) MatZeroEntries(H);
199: /* Get a pointer to vector data */
200: VecGetArrayRead(X,&x);
202: /* Compute H(X) entries */
203: if (user->chained) {
204: for (i=0; i<user->n-1; i++) {
205: PetscScalar t1 = x[i+1] - x[i]*x[i];
206: v[0][0] = 2 + 2*alpha*(t1*(-2) - 2*x[i]);
207: v[0][1] = 2*alpha*(-2*x[i]);
208: v[1][0] = 2*alpha*(-2*x[i]);
209: v[1][1] = 2*alpha*t1;
210: ind[0] = i; ind[1] = i+1;
211: MatSetValues(H,2,ind,2,ind,v[0],ADD_VALUES);
212: }
213: } else {
214: for (i=0; i<user->n/2; i++) {
215: v[1][1] = 2*alpha;
216: v[0][0] = -4*alpha*(x[2*i+1]-3*x[2*i]*x[2*i]) + 2;
217: v[1][0] = v[0][1] = -4.0*alpha*x[2*i];
218: ind[0]=2*i; ind[1]=2*i+1;
219: MatSetValues(H,2,ind,2,ind,v[0],INSERT_VALUES);
220: }
221: }
222: VecRestoreArrayRead(X,&x);
224: /* Assemble matrix */
225: MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
226: MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
227: PetscLogFlops(9.0*user->n/2.0);
228: return 0;
229: }
231: /*TEST
233: build:
234: requires: !complex
236: test:
237: args: -tao_type bqnls -tao_monitor
238: requires: !single
240: TEST*/