1 /* Program usage: mpiexec -n 1 rosenbrock1 [-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 KSP ksp;
38 PC pc;
39 Mat M;
40 Vec in, out, out2;
41 PetscReal mult_solve_dist;
42
43 /* Initialize TAO and PETSc */
44 PetscFunctionBeginUser;
45 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
46 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
47 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Incorrect number of processors");
48
49 /* Initialize problem parameters */
50 user.n = 2;
51 user.alpha = 99.0;
52 user.chained = PETSC_FALSE;
53 /* Check for command line arguments to override defaults */
54 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &user.n, &flg));
55 PetscCall(PetscOptionsGetReal(NULL, NULL, "-alpha", &user.alpha, &flg));
56 PetscCall(PetscOptionsGetBool(NULL, NULL, "-chained", &user.chained, &flg));
57 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_lmvm", &test_lmvm, &flg));
58
59 /* Allocate vectors for the solution and gradient */
60 PetscCall(VecCreateSeq(PETSC_COMM_SELF, user.n, &x));
61 PetscCall(MatCreateSeqBAIJ(PETSC_COMM_SELF, 2, user.n, user.n, 1, NULL, &H));
62
63 /* The TAO code begins here */
64
65 /* Create TAO solver with desired solution method */
66 PetscCall(TaoCreate(PETSC_COMM_SELF, &tao));
67 PetscCall(TaoSetType(tao, TAOLMVM));
68
69 /* Set solution vec and an initial guess */
70 PetscCall(VecSet(x, zero));
71 PetscCall(TaoSetSolution(tao, x));
72
73 /* Set routines for function, gradient, hessian evaluation */
74 PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, &user));
75 PetscCall(TaoSetHessian(tao, H, H, FormHessian, &user));
76
77 /* Test the LMVM matrix */
78 if (test_lmvm) PetscCall(PetscOptionsSetValue(NULL, "-tao_type", "bqnktr"));
79
80 /* Check for TAO command line options */
81 PetscCall(TaoSetFromOptions(tao));
82
83 /* SOLVE THE APPLICATION */
84 PetscCall(TaoSolve(tao));
85
86 /* Test the LMVM matrix */
87 if (test_lmvm) {
88 PetscCall(TaoGetKSP(tao, &ksp));
89 PetscCall(KSPGetPC(ksp, &pc));
90 PetscCall(PCLMVMGetMatLMVM(pc, &M));
91 PetscCall(VecDuplicate(x, &in));
92 PetscCall(VecDuplicate(x, &out));
93 PetscCall(VecDuplicate(x, &out2));
94 PetscCall(VecSet(in, 1.0));
95 PetscCall(MatMult(M, in, out));
96 PetscCall(MatSolve(M, out, out2));
97 PetscCall(VecAXPY(out2, -1.0, in));
98 PetscCall(VecNorm(out2, NORM_2, &mult_solve_dist));
99 if (mult_solve_dist < 1.e-11) {
100 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)tao), "Inverse error of LMVM MatMult and MatSolve: < 1.e-11\n"));
101 } else if (mult_solve_dist < 1.e-6) {
102 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)tao), "Inverse error of LMVM MatMult and MatSolve: < 1.e-6\n"));
103 } else {
104 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)tao), "Inverse error of LMVM MatMult and MatSolve is not small: %e\n", (double)mult_solve_dist));
105 }
106 PetscCall(VecDestroy(&in));
107 PetscCall(VecDestroy(&out));
108 PetscCall(VecDestroy(&out2));
109 }
110
111 PetscCall(TaoDestroy(&tao));
112 PetscCall(VecDestroy(&x));
113 PetscCall(MatDestroy(&H));
114
115 PetscCall(PetscFinalize());
116 return 0;
117 }
118
119 /* -------------------------------------------------------------------- */
120 /*
121 FormFunctionGradient - Evaluates the function, f(X), and gradient, G(X).
122
123 Input Parameters:
124 . tao - the Tao context
125 . X - input vector
126 . ptr - optional user-defined context, as set by TaoSetFunctionGradient()
127
128 Output Parameters:
129 . G - vector containing the newly evaluated gradient
130 . f - function value
131
132 Note:
133 Some optimization methods ask for the function and the gradient evaluation
134 at the same time. Evaluating both at once may be more efficient that
135 evaluating each separately.
136 */
FormFunctionGradient(Tao tao,Vec X,PetscReal * f,Vec G,void * ptr)137 PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
138 {
139 AppCtx *user = (AppCtx *)ptr;
140 PetscInt i, nn = user->n / 2;
141 PetscReal ff = 0, t1, t2, alpha = user->alpha;
142 PetscScalar *g;
143 const PetscScalar *x;
144
145 PetscFunctionBeginUser;
146 /* Get pointers to vector data */
147 PetscCall(VecGetArrayRead(X, &x));
148 PetscCall(VecGetArray(G, &g));
149
150 /* Compute G(X) */
151 if (user->chained) {
152 g[0] = 0;
153 for (i = 0; i < user->n - 1; i++) {
154 t1 = x[i + 1] - x[i] * x[i];
155 ff += PetscSqr(1 - x[i]) + alpha * t1 * t1;
156 g[i] += -2 * (1 - x[i]) + 2 * alpha * t1 * (-2 * x[i]);
157 g[i + 1] = 2 * alpha * t1;
158 }
159 } else {
160 for (i = 0; i < nn; i++) {
161 t1 = x[2 * i + 1] - x[2 * i] * x[2 * i];
162 t2 = 1 - x[2 * i];
163 ff += alpha * t1 * t1 + t2 * t2;
164 g[2 * i] = -4 * alpha * t1 * x[2 * i] - 2.0 * t2;
165 g[2 * i + 1] = 2 * alpha * t1;
166 }
167 }
168
169 /* Restore vectors */
170 PetscCall(VecRestoreArrayRead(X, &x));
171 PetscCall(VecRestoreArray(G, &g));
172 *f = ff;
173
174 PetscCall(PetscLogFlops(15.0 * nn));
175 PetscFunctionReturn(PETSC_SUCCESS);
176 }
177
178 /* ------------------------------------------------------------------- */
179 /*
180 FormHessian - Evaluates Hessian matrix.
181
182 Input Parameters:
183 . tao - the Tao context
184 . x - input vector
185 . ptr - optional user-defined context, as set by TaoSetHessian()
186
187 Output Parameters:
188 . H - Hessian matrix
189
190 Note: Providing the Hessian may not be necessary. Only some solvers
191 require this matrix.
192 */
FormHessian(Tao tao,Vec X,Mat H,Mat Hpre,void * ptr)193 PetscErrorCode FormHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr)
194 {
195 AppCtx *user = (AppCtx *)ptr;
196 PetscInt i, ind[2];
197 PetscReal alpha = user->alpha;
198 PetscReal v[2][2];
199 const PetscScalar *x;
200 PetscBool assembled;
201
202 PetscFunctionBeginUser;
203 /* Zero existing matrix entries */
204 PetscCall(MatAssembled(H, &assembled));
205 if (assembled) PetscCall(MatZeroEntries(H));
206
207 /* Get a pointer to vector data */
208 PetscCall(VecGetArrayRead(X, &x));
209
210 /* Compute H(X) entries */
211 if (user->chained) {
212 PetscCall(MatZeroEntries(H));
213 for (i = 0; i < user->n - 1; i++) {
214 PetscScalar t1 = x[i + 1] - x[i] * x[i];
215 v[0][0] = 2 + 2 * alpha * (t1 * (-2) - 2 * x[i]);
216 v[0][1] = 2 * alpha * (-2 * x[i]);
217 v[1][0] = 2 * alpha * (-2 * x[i]);
218 v[1][1] = 2 * alpha * t1;
219 ind[0] = i;
220 ind[1] = i + 1;
221 PetscCall(MatSetValues(H, 2, ind, 2, ind, v[0], ADD_VALUES));
222 }
223 } else {
224 for (i = 0; i < user->n / 2; i++) {
225 v[1][1] = 2 * alpha;
226 v[0][0] = -4 * alpha * (x[2 * i + 1] - 3 * x[2 * i] * x[2 * i]) + 2;
227 v[1][0] = v[0][1] = -4.0 * alpha * x[2 * i];
228 ind[0] = 2 * i;
229 ind[1] = 2 * i + 1;
230 PetscCall(MatSetValues(H, 2, ind, 2, ind, v[0], INSERT_VALUES));
231 }
232 }
233 PetscCall(VecRestoreArrayRead(X, &x));
234
235 /* Assemble matrix */
236 PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
237 PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
238 PetscCall(PetscLogFlops(9.0 * user->n / 2.0));
239 PetscFunctionReturn(PETSC_SUCCESS);
240 }
241
242 /*TEST
243
244 build:
245 requires: !complex
246
247 test:
248 requires: !single
249 args: -tao_monitor_short -tao_type nls -tao_gatol 1.e-4
250
251 test:
252 suffix: 2
253 requires: !single
254 args: -tao_monitor_short -tao_type lmvm -tao_gatol 1.e-3
255
256 test:
257 suffix: 3
258 requires: !single
259 args: -tao_monitor_short -tao_type ntr -tao_gatol 1.e-4
260
261 test:
262 suffix: 4
263 requires: !single
264 args: -tao_monitor_short -tao_type ntr -tao_mf_hessian -tao_ntr_pc_type none -tao_gatol 1.e-4
265
266 test:
267 suffix: 5
268 requires: !single
269 args: -tao_monitor_short -tao_type bntr -tao_gatol 1.e-4
270
271 test:
272 suffix: 6
273 requires: !single
274 args: -tao_monitor_short -tao_type bntl -tao_gatol 1.e-4
275
276 test:
277 suffix: 7
278 requires: !single
279 args: -tao_monitor_short -tao_type bnls -tao_gatol 1.e-4
280
281 test:
282 suffix: 8
283 requires: !single
284 args: -tao_monitor_short -tao_type bntr -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 -tao_bnk_cg_tao_monitor_short
285
286 test:
287 suffix: 9
288 requires: !single
289 args: -tao_monitor_short -tao_type bntl -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 -tao_bnk_cg_tao_monitor_short
290
291 test:
292 suffix: 10
293 requires: !single
294 args: -tao_monitor_short -tao_type bnls -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 -tao_bnk_cg_tao_monitor_short
295
296 test:
297 suffix: 11
298 requires: !single
299 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmbroyden
300
301 test:
302 suffix: 12
303 requires: !single
304 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmbadbroyden
305
306 test:
307 suffix: 13
308 requires: !single
309 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmsymbroyden -tao_bqnk_mat_lmvm_beta {{0.0 0.25 1.0}} -tao_bqnk_mat_lmvm_rho 0.75 -tao_bqnk_mat_lmvm_sigma_hist 2
310
311 test:
312 suffix: 14
313 requires: !single
314 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmbfgs -tao_bqnk_mat_lmvm_scale_type {{scalar diagonal}} -tao_bqnk_mat_lmvm_alpha {{0.0 0.25 0.5}} -tao_bqnk_mat_lmvm_theta 1.0
315
316 test:
317 suffix: 15
318 requires: !single
319 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmdfp
320
321 test:
322 suffix: 16
323 requires: !single
324 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmsr1
325
326 test:
327 suffix: 17
328 requires: !single
329 args: -tao_monitor_short -tao_gatol 1e-4 -tao_type bqnls
330
331 test:
332 suffix: 18
333 requires: !single
334 args: -tao_monitor_short -tao_gatol 1e-4 -tao_type blmvm
335
336 test:
337 suffix: 19
338 requires: !single
339 args: -tao_monitor_short -tao_gatol 1e-4 -tao_type bqnktr -tao_bqnk_mat_type lmvmsr1
340
341 test:
342 suffix: 20
343 requires: !single
344 args: -tao_monitor -tao_gatol 1e-4 -tao_type blmvm -tao_ls_monitor
345
346 test:
347 suffix: 21
348 requires: !single
349 args: -test_lmvm -tao_max_it 10 -tao_bqnk_mat_type lmvmsymbadbroyden
350
351 test:
352 suffix: 22
353 requires: !single
354 args: -tao_max_it 1 -tao_converged_reason
355
356 test:
357 suffix: 23
358 requires: !single
359 args: -tao_max_funcs 0 -tao_converged_reason
360
361 test:
362 suffix: 24
363 requires: !single
364 args: -tao_gatol 10 -tao_converged_reason
365
366 test:
367 suffix: 25
368 requires: !single
369 args: -tao_grtol 10 -tao_converged_reason
370
371 test:
372 suffix: 26
373 requires: !single
374 args: -tao_gttol 10 -tao_converged_reason
375
376 test:
377 suffix: 27
378 requires: !single
379 args: -tao_steptol 10 -tao_converged_reason
380
381 test:
382 suffix: 28
383 requires: !single
384 args: -tao_fmin 10 -tao_converged_reason
385
386 test:
387 suffix: snes
388 requires: !single
389 args: -snes_monitor ::ascii_info_detail -tao_type snes -snes_type newtontr -snes_atol 1.e-4 -pc_type none -tao_mf_hessian -ksp_type cg
390
391 test:
392 suffix: snes_ls_armijo
393 requires: !single
394 args: -snes_monitor ::ascii_info_detail -tao_type snes -snes_type newtonls -snes_atol 1.e-4 -pc_type none -tao_mf_hessian -snes_linesearch_monitor -snes_linesearch_order 1
395
396 test:
397 suffix: snes_tr_cgnegcurve_kmdc
398 requires: !single
399 args: -snes_monitor ::ascii_info_detail -tao_type snes -snes_type newtontr -snes_atol 1.e-4 -pc_type none -ksp_type cg -snes_tr_kmdc 0.9 -ksp_converged_neg_curve -ksp_converged_reason
400
401 test:
402 suffix: snes_ls_lmvm
403 requires: !single
404 args: -snes_monitor ::ascii_info_detail -tao_type snes -snes_type newtonls -snes_atol 1.e-4 -pc_type lmvm -tao_mf_hessian
405
406 TEST*/
407