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, TAOLMVM));
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(TaoLMVMRecycle(tao, PETSC_TRUE));
81 reason = TAO_CONTINUE_ITERATING;
82 while (reason != TAO_CONVERGED_GATOL) {
83 PetscCall(TaoSolve(tao));
84 PetscCall(TaoGetConvergedReason(tao, &reason));
85 PetscCall(TaoGetIterationNumber(tao, &its));
86 recycled_its += its;
87 PetscCall(PetscPrintf(PETSC_COMM_SELF, "-----------------------\n"));
88 }
89
90 /* Disable recycling and solve again! */
91 PetscCall(TaoSetMaximumIterations(tao, 100));
92 PetscCall(TaoLMVMRecycle(tao, PETSC_FALSE));
93 PetscCall(VecSet(x, zero));
94 PetscCall(TaoSolve(tao));
95 PetscCall(TaoGetConvergedReason(tao, &reason));
96 PetscCheck(reason == TAO_CONVERGED_GATOL, PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "Solution failed to converge!");
97 PetscCall(TaoGetIterationNumber(tao, &oneshot_its));
98 PetscCall(PetscPrintf(PETSC_COMM_SELF, "-----------------------\n"));
99 PetscCall(PetscPrintf(PETSC_COMM_SELF, "recycled its: %" PetscInt_FMT " | oneshot its: %" PetscInt_FMT "\n", recycled_its, oneshot_its));
100 PetscCheck(recycled_its == oneshot_its, PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "LMVM recycling does not work!");
101
102 PetscCall(TaoDestroy(&tao));
103 PetscCall(VecDestroy(&x));
104 PetscCall(MatDestroy(&H));
105
106 PetscCall(PetscFinalize());
107 return 0;
108 }
109
110 /* -------------------------------------------------------------------- */
111 /*
112 FormFunctionGradient - Evaluates the function, f(X), and gradient, G(X).
113
114 Input Parameters:
115 . tao - the Tao context
116 . X - input vector
117 . ptr - optional user-defined context, as set by TaoSetFunctionGradient()
118
119 Output Parameters:
120 . G - vector containing the newly evaluated gradient
121 . f - function value
122
123 Note:
124 Some optimization methods ask for the function and the gradient evaluation
125 at the same time. Evaluating both at once may be more efficient that
126 evaluating each separately.
127 */
FormFunctionGradient(Tao tao,Vec X,PetscReal * f,Vec G,void * ptr)128 PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
129 {
130 AppCtx *user = (AppCtx *)ptr;
131 PetscInt i, nn = user->n / 2;
132 PetscReal ff = 0, t1, t2, alpha = user->alpha;
133 PetscScalar *g;
134 const PetscScalar *x;
135
136 PetscFunctionBeginUser;
137 /* Get pointers to vector data */
138 PetscCall(VecGetArrayRead(X, &x));
139 PetscCall(VecGetArray(G, &g));
140
141 /* Compute G(X) */
142 if (user->chained) {
143 g[0] = 0;
144 for (i = 0; i < user->n - 1; i++) {
145 t1 = x[i + 1] - x[i] * x[i];
146 ff += PetscSqr(1 - x[i]) + alpha * t1 * t1;
147 g[i] += -2 * (1 - x[i]) + 2 * alpha * t1 * (-2 * x[i]);
148 g[i + 1] = 2 * alpha * t1;
149 }
150 } else {
151 for (i = 0; i < nn; i++) {
152 t1 = x[2 * i + 1] - x[2 * i] * x[2 * i];
153 t2 = 1 - x[2 * i];
154 ff += alpha * t1 * t1 + t2 * t2;
155 g[2 * i] = -4 * alpha * t1 * x[2 * i] - 2.0 * t2;
156 g[2 * i + 1] = 2 * alpha * t1;
157 }
158 }
159
160 /* Restore vectors */
161 PetscCall(VecRestoreArrayRead(X, &x));
162 PetscCall(VecRestoreArray(G, &g));
163 *f = ff;
164
165 PetscCall(PetscLogFlops(15.0 * nn));
166 PetscFunctionReturn(PETSC_SUCCESS);
167 }
168
169 /* ------------------------------------------------------------------- */
170 /*
171 FormHessian - Evaluates Hessian matrix.
172
173 Input Parameters:
174 . tao - the Tao context
175 . x - input vector
176 . ptr - optional user-defined context, as set by TaoSetHessian()
177
178 Output Parameters:
179 . H - Hessian matrix
180
181 Note: Providing the Hessian may not be necessary. Only some solvers
182 require this matrix.
183 */
FormHessian(Tao tao,Vec X,Mat H,Mat Hpre,void * ptr)184 PetscErrorCode FormHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr)
185 {
186 AppCtx *user = (AppCtx *)ptr;
187 PetscInt i, ind[2];
188 PetscReal alpha = user->alpha;
189 PetscReal v[2][2];
190 const PetscScalar *x;
191 PetscBool assembled;
192
193 PetscFunctionBeginUser;
194 /* Zero existing matrix entries */
195 PetscCall(MatAssembled(H, &assembled));
196 if (assembled) PetscCall(MatZeroEntries(H));
197
198 /* Get a pointer to vector data */
199 PetscCall(VecGetArrayRead(X, &x));
200
201 /* Compute H(X) entries */
202 if (user->chained) {
203 PetscCall(MatZeroEntries(H));
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;
211 ind[1] = i + 1;
212 PetscCall(MatSetValues(H, 2, ind, 2, ind, v[0], ADD_VALUES));
213 }
214 } else {
215 for (i = 0; i < user->n / 2; i++) {
216 v[1][1] = 2 * alpha;
217 v[0][0] = -4 * alpha * (x[2 * i + 1] - 3 * x[2 * i] * x[2 * i]) + 2;
218 v[1][0] = v[0][1] = -4.0 * alpha * x[2 * i];
219 ind[0] = 2 * i;
220 ind[1] = 2 * i + 1;
221 PetscCall(MatSetValues(H, 2, ind, 2, ind, v[0], INSERT_VALUES));
222 }
223 }
224 PetscCall(VecRestoreArrayRead(X, &x));
225
226 /* Assemble matrix */
227 PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
228 PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
229 PetscCall(PetscLogFlops(9.0 * user->n / 2.0));
230 PetscFunctionReturn(PETSC_SUCCESS);
231 }
232
233 /*TEST
234
235 build:
236 requires: !complex
237
238 test:
239 args: -tao_type lmvm -tao_monitor
240 requires: !single
241
242 TEST*/
243