1 /* Program usage: mpiexec -n 1 rosenbrock2 [-help] [all TAO options] */
2
3 /* Include "petsctao.h" so we can use TAO solvers. */
4 #include <petsctao.h>
5
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";
12
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;
23
24 /* -------------- User-defined routines ---------- */
25 PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
26 PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *);
27
main(int argc,char ** argv)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;
39
40 /* Initialize TAO and PETSc */
41 PetscFunctionBeginUser;
42 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
43 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
44 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Incorrect number of processors");
45
46 /* Initialize problem parameters */
47 user.n = 2;
48 user.alpha = 99.0;
49 user.chained = PETSC_FALSE;
50 /* Check for command line arguments to override defaults */
51 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &user.n, &flg));
52 PetscCall(PetscOptionsGetReal(NULL, NULL, "-alpha", &user.alpha, &flg));
53 PetscCall(PetscOptionsGetBool(NULL, NULL, "-chained", &user.chained, &flg));
54 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_lmvm", &test_lmvm, &flg));
55
56 /* Allocate vectors for the solution and gradient */
57 PetscCall(VecCreateSeq(PETSC_COMM_SELF, user.n, &x));
58 PetscCall(MatCreateSeqBAIJ(PETSC_COMM_SELF, 2, user.n, user.n, 1, NULL, &H));
59
60 /* The TAO code begins here */
61
62 /* Create TAO solver with desired solution method */
63 PetscCall(TaoCreate(PETSC_COMM_SELF, &tao));
64 PetscCall(TaoSetType(tao, TAOBQNLS));
65
66 /* Set solution vec and an initial guess */
67 PetscCall(VecSet(x, zero));
68 PetscCall(TaoSetSolution(tao, x));
69
70 /* Set routines for function, gradient, hessian evaluation */
71 PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, &user));
72 PetscCall(TaoSetHessian(tao, H, H, FormHessian, &user));
73
74 /* Check for TAO command line options */
75 PetscCall(TaoSetFromOptions(tao));
76
77 /* Solve the problem */
78 PetscCall(TaoSetTolerances(tao, 1.e-5, 0.0, 0.0));
79 PetscCall(TaoSetMaximumIterations(tao, 5));
80 PetscCall(TaoSetRecycleHistory(tao, PETSC_TRUE));
81 reason = TAO_CONTINUE_ITERATING;
82 flg = PETSC_FALSE;
83 PetscCall(TaoGetRecycleHistory(tao, &flg));
84 if (flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Recycle: enabled\n"));
85 while (reason != TAO_CONVERGED_GATOL) {
86 PetscCall(TaoSolve(tao));
87 PetscCall(TaoGetConvergedReason(tao, &reason));
88 PetscCall(TaoGetIterationNumber(tao, &its));
89 recycled_its += its;
90 PetscCall(PetscPrintf(PETSC_COMM_SELF, "-----------------------\n"));
91 }
92
93 /* Disable recycling and solve again! */
94 PetscCall(TaoSetMaximumIterations(tao, 100));
95 PetscCall(TaoSetRecycleHistory(tao, PETSC_FALSE));
96 PetscCall(VecSet(x, zero));
97 PetscCall(TaoGetRecycleHistory(tao, &flg));
98 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Recycle: disabled\n"));
99 PetscCall(TaoSolve(tao));
100 PetscCall(TaoGetConvergedReason(tao, &reason));
101 PetscCheck(reason == TAO_CONVERGED_GATOL, PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "Solution failed to converge!");
102 PetscCall(TaoGetIterationNumber(tao, &oneshot_its));
103 PetscCall(PetscPrintf(PETSC_COMM_SELF, "-----------------------\n"));
104 PetscCall(PetscPrintf(PETSC_COMM_SELF, "recycled its: %" PetscInt_FMT " | oneshot its: %" PetscInt_FMT "\n", recycled_its, oneshot_its));
105 PetscCheck(recycled_its == oneshot_its, PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "Recycled solution does not match oneshot solution!");
106
107 PetscCall(TaoDestroy(&tao));
108 PetscCall(VecDestroy(&x));
109 PetscCall(MatDestroy(&H));
110
111 PetscCall(PetscFinalize());
112 return 0;
113 }
114
115 /* -------------------------------------------------------------------- */
116 /*
117 FormFunctionGradient - Evaluates the function, f(X), and gradient, G(X).
118
119 Input Parameters:
120 . tao - the Tao context
121 . X - input vector
122 . ptr - optional user-defined context, as set by TaoSetFunctionGradient()
123
124 Output Parameters:
125 . G - vector containing the newly evaluated gradient
126 . f - function value
127
128 Note:
129 Some optimization methods ask for the function and the gradient evaluation
130 at the same time. Evaluating both at once may be more efficient than
131 evaluating each separately.
132 */
FormFunctionGradient(Tao tao,Vec X,PetscReal * f,Vec G,void * ptr)133 PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
134 {
135 AppCtx *user = (AppCtx *)ptr;
136 PetscInt i, nn = user->n / 2;
137 PetscReal ff = 0, t1, t2, alpha = user->alpha;
138 PetscScalar *g;
139 const PetscScalar *x;
140
141 PetscFunctionBeginUser;
142 /* Get pointers to vector data */
143 PetscCall(VecGetArrayRead(X, &x));
144 PetscCall(VecGetArrayWrite(G, &g));
145
146 /* Compute G(X) */
147 if (user->chained) {
148 g[0] = 0;
149 for (i = 0; i < user->n - 1; i++) {
150 t1 = x[i + 1] - x[i] * x[i];
151 ff += PetscSqr(1 - x[i]) + alpha * t1 * t1;
152 g[i] += -2 * (1 - x[i]) + 2 * alpha * t1 * (-2 * x[i]);
153 g[i + 1] = 2 * alpha * t1;
154 }
155 } else {
156 for (i = 0; i < nn; i++) {
157 t1 = x[2 * i + 1] - x[2 * i] * x[2 * i];
158 t2 = 1 - x[2 * i];
159 ff += alpha * t1 * t1 + t2 * t2;
160 g[2 * i] = -4 * alpha * t1 * x[2 * i] - 2.0 * t2;
161 g[2 * i + 1] = 2 * alpha * t1;
162 }
163 }
164
165 /* Restore vectors */
166 PetscCall(VecRestoreArrayRead(X, &x));
167 PetscCall(VecRestoreArrayWrite(G, &g));
168 *f = ff;
169
170 PetscCall(PetscLogFlops(15.0 * nn));
171 PetscFunctionReturn(PETSC_SUCCESS);
172 }
173
174 /* ------------------------------------------------------------------- */
175 /*
176 FormHessian - Evaluates Hessian matrix.
177
178 Input Parameters:
179 . tao - the Tao context
180 . x - input vector
181 . ptr - optional user-defined context, as set by TaoSetHessian()
182
183 Output Parameters:
184 . H - Hessian matrix
185
186 Note: Providing the Hessian may not be necessary. Only some solvers
187 require this matrix.
188 */
FormHessian(Tao tao,Vec X,Mat H,Mat Hpre,void * ptr)189 PetscErrorCode FormHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr)
190 {
191 AppCtx *user = (AppCtx *)ptr;
192 PetscInt i, ind[2];
193 PetscReal alpha = user->alpha;
194 PetscReal v[2][2];
195 const PetscScalar *x;
196 PetscBool assembled;
197
198 PetscFunctionBeginUser;
199 /* Zero existing matrix entries */
200 PetscCall(MatAssembled(H, &assembled));
201 if (assembled || user->chained) PetscCall(MatZeroEntries(H));
202
203 /* Get a pointer to vector data */
204 PetscCall(VecGetArrayRead(X, &x));
205
206 /* Compute H(X) entries */
207 if (user->chained) {
208 for (i = 0; i < user->n - 1; i++) {
209 PetscScalar t1 = x[i + 1] - x[i] * x[i];
210 v[0][0] = 2 + 2 * alpha * (t1 * (-2) - 2 * x[i]);
211 v[0][1] = 2 * alpha * (-2 * x[i]);
212 v[1][0] = 2 * alpha * (-2 * x[i]);
213 v[1][1] = 2 * alpha * t1;
214 ind[0] = i;
215 ind[1] = i + 1;
216 PetscCall(MatSetValues(H, 2, ind, 2, ind, v[0], ADD_VALUES));
217 }
218 } else {
219 for (i = 0; i < user->n / 2; i++) {
220 v[1][1] = 2 * alpha;
221 v[0][0] = -4 * alpha * (x[2 * i + 1] - 3 * x[2 * i] * x[2 * i]) + 2;
222 v[1][0] = v[0][1] = -4.0 * alpha * x[2 * i];
223 ind[0] = 2 * i;
224 ind[1] = 2 * i + 1;
225 PetscCall(MatSetValues(H, 2, ind, 2, ind, v[0], INSERT_VALUES));
226 }
227 }
228 PetscCall(VecRestoreArrayRead(X, &x));
229
230 /* Assemble matrix */
231 PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
232 PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
233 PetscCall(PetscLogFlops(9.0 * user->n / 2.0));
234 PetscFunctionReturn(PETSC_SUCCESS);
235 }
236
237 /*TEST
238
239 build:
240 requires: !complex
241
242 test:
243 args: -tao_type bqnls -tao_monitor
244 requires: !single
245
246 TEST*/
247