1 #include <../src/tao/complementarity/impls/ssls/ssls.h>
2 /*
3 Context for ASXLS
4 -- active-set - reduced matrices formed
5 - inherit properties of original system
6 -- semismooth (S) - function not differentiable
7 - merit function continuously differentiable
8 - Fischer-Burmeister reformulation of complementarity
9 - Billups composition for two finite bounds
10 -- infeasible (I) - iterates not guaranteed to remain within bounds
11 -- feasible (F) - iterates guaranteed to remain within bounds
12 -- linesearch (LS) - Armijo rule on direction
13
14 Many other reformulations are possible and combinations of
15 feasible/infeasible and linesearch/trust region are possible.
16
17 Basic theory
18 Fischer-Burmeister reformulation is semismooth with a continuously
19 differentiable merit function and strongly semismooth if the F has
20 lipschitz continuous derivatives.
21
22 Every accumulation point generated by the algorithm is a stationary
23 point for the merit function. Stationary points of the merit function
24 are solutions of the complementarity problem if
25 a. the stationary point has a BD-regular subdifferential, or
26 b. the Schur complement F'/F'_ff is a P_0-matrix where ff is the
27 index set corresponding to the free variables.
28
29 If one of the accumulation points has a BD-regular subdifferential then
30 a. the entire sequence converges to this accumulation point at
31 a local q-superlinear rate
32 b. if in addition the reformulation is strongly semismooth near
33 this accumulation point, then the algorithm converges at a
34 local q-quadratic rate.
35
36 The theory for the feasible version follows from the feasible descent
37 algorithm framework.
38
39 References:
40 Billups, "Algorithms for Complementarity Problems and Generalized
41 Equations," Ph.D thesis, University of Wisconsin Madison, 1995.
42 De Luca, Facchinei, Kanzow, "A Semismooth Equation Approach to the
43 Solution of Nonlinear Complementarity Problems," Mathematical
44 Programming, 75, pages 407439, 1996.
45 Ferris, Kanzow, Munson, "Feasible Descent Algorithms for Mixed
46 Complementarity Problems," Mathematical Programming, 86,
47 pages 475497, 1999.
48 Fischer, "A Special Newton type Optimization Method," Optimization,
49 24, 1992
50 Munson, Facchinei, Ferris, Fischer, Kanzow, "The Semismooth Algorithm
51 for Large Scale Complementarity Problems," Technical Report,
52 University of Wisconsin Madison, 1999.
53 */
54
55
TaoSetUp_ASFLS(Tao tao)56 static PetscErrorCode TaoSetUp_ASFLS(Tao tao)
57 {
58 TAO_SSLS *asls = (TAO_SSLS *)tao->data;
59 PetscErrorCode ierr;
60
61 PetscFunctionBegin;
62 ierr = VecDuplicate(tao->solution,&tao->gradient);CHKERRQ(ierr);
63 ierr = VecDuplicate(tao->solution,&tao->stepdirection);CHKERRQ(ierr);
64 ierr = VecDuplicate(tao->solution,&asls->ff);CHKERRQ(ierr);
65 ierr = VecDuplicate(tao->solution,&asls->dpsi);CHKERRQ(ierr);
66 ierr = VecDuplicate(tao->solution,&asls->da);CHKERRQ(ierr);
67 ierr = VecDuplicate(tao->solution,&asls->db);CHKERRQ(ierr);
68 ierr = VecDuplicate(tao->solution,&asls->t1);CHKERRQ(ierr);
69 ierr = VecDuplicate(tao->solution,&asls->t2);CHKERRQ(ierr);
70 ierr = VecDuplicate(tao->solution, &asls->w);CHKERRQ(ierr);
71 asls->fixed = NULL;
72 asls->free = NULL;
73 asls->J_sub = NULL;
74 asls->Jpre_sub = NULL;
75 asls->r1 = NULL;
76 asls->r2 = NULL;
77 asls->r3 = NULL;
78 asls->dxfree = NULL;
79 PetscFunctionReturn(0);
80 }
81
Tao_ASLS_FunctionGradient(TaoLineSearch ls,Vec X,PetscReal * fcn,Vec G,void * ptr)82 static PetscErrorCode Tao_ASLS_FunctionGradient(TaoLineSearch ls, Vec X, PetscReal *fcn, Vec G, void *ptr)
83 {
84 Tao tao = (Tao)ptr;
85 TAO_SSLS *asls = (TAO_SSLS *)tao->data;
86 PetscErrorCode ierr;
87
88 PetscFunctionBegin;
89 ierr = TaoComputeConstraints(tao, X, tao->constraints);CHKERRQ(ierr);
90 ierr = VecFischer(X,tao->constraints,tao->XL,tao->XU,asls->ff);CHKERRQ(ierr);
91 ierr = VecNorm(asls->ff,NORM_2,&asls->merit);CHKERRQ(ierr);
92 *fcn = 0.5*asls->merit*asls->merit;
93 ierr = TaoComputeJacobian(tao,tao->solution,tao->jacobian,tao->jacobian_pre);CHKERRQ(ierr);
94
95 ierr = MatDFischer(tao->jacobian, tao->solution, tao->constraints,tao->XL, tao->XU, asls->t1, asls->t2,asls->da, asls->db);CHKERRQ(ierr);
96 ierr = VecPointwiseMult(asls->t1, asls->ff, asls->db);CHKERRQ(ierr);
97 ierr = MatMultTranspose(tao->jacobian,asls->t1,G);CHKERRQ(ierr);
98 ierr = VecPointwiseMult(asls->t1, asls->ff, asls->da);CHKERRQ(ierr);
99 ierr = VecAXPY(G,1.0,asls->t1);CHKERRQ(ierr);
100 PetscFunctionReturn(0);
101 }
102
TaoDestroy_ASFLS(Tao tao)103 static PetscErrorCode TaoDestroy_ASFLS(Tao tao)
104 {
105 TAO_SSLS *ssls = (TAO_SSLS *)tao->data;
106 PetscErrorCode ierr;
107
108 PetscFunctionBegin;
109 ierr = VecDestroy(&ssls->ff);CHKERRQ(ierr);
110 ierr = VecDestroy(&ssls->dpsi);CHKERRQ(ierr);
111 ierr = VecDestroy(&ssls->da);CHKERRQ(ierr);
112 ierr = VecDestroy(&ssls->db);CHKERRQ(ierr);
113 ierr = VecDestroy(&ssls->w);CHKERRQ(ierr);
114 ierr = VecDestroy(&ssls->t1);CHKERRQ(ierr);
115 ierr = VecDestroy(&ssls->t2);CHKERRQ(ierr);
116 ierr = VecDestroy(&ssls->r1);CHKERRQ(ierr);
117 ierr = VecDestroy(&ssls->r2);CHKERRQ(ierr);
118 ierr = VecDestroy(&ssls->r3);CHKERRQ(ierr);
119 ierr = VecDestroy(&ssls->dxfree);CHKERRQ(ierr);
120 ierr = MatDestroy(&ssls->J_sub);CHKERRQ(ierr);
121 ierr = MatDestroy(&ssls->Jpre_sub);CHKERRQ(ierr);
122 ierr = ISDestroy(&ssls->fixed);CHKERRQ(ierr);
123 ierr = ISDestroy(&ssls->free);CHKERRQ(ierr);
124 ierr = PetscFree(tao->data);CHKERRQ(ierr);
125 tao->data = NULL;
126 PetscFunctionReturn(0);
127 }
128
TaoSolve_ASFLS(Tao tao)129 static PetscErrorCode TaoSolve_ASFLS(Tao tao)
130 {
131 TAO_SSLS *asls = (TAO_SSLS *)tao->data;
132 PetscReal psi,ndpsi, normd, innerd, t=0;
133 PetscInt nf;
134 PetscErrorCode ierr;
135 TaoLineSearchConvergedReason ls_reason;
136
137 PetscFunctionBegin;
138 /* Assume that Setup has been called!
139 Set the structure for the Jacobian and create a linear solver. */
140
141 ierr = TaoComputeVariableBounds(tao);CHKERRQ(ierr);
142 ierr = TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,Tao_ASLS_FunctionGradient,tao);CHKERRQ(ierr);
143 ierr = TaoLineSearchSetObjectiveRoutine(tao->linesearch,Tao_SSLS_Function,tao);CHKERRQ(ierr);
144 ierr = TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);CHKERRQ(ierr);
145
146 ierr = VecMedian(tao->XL, tao->solution, tao->XU, tao->solution);CHKERRQ(ierr);
147
148 /* Calculate the function value and fischer function value at the
149 current iterate */
150 ierr = TaoLineSearchComputeObjectiveAndGradient(tao->linesearch,tao->solution,&psi,asls->dpsi);CHKERRQ(ierr);
151 ierr = VecNorm(asls->dpsi,NORM_2,&ndpsi);CHKERRQ(ierr);
152
153 tao->reason = TAO_CONTINUE_ITERATING;
154 while (1) {
155 /* Check the converged criteria */
156 ierr = PetscInfo3(tao,"iter %D, merit: %g, ||dpsi||: %g\n",tao->niter,(double)asls->merit,(double)ndpsi);CHKERRQ(ierr);
157 ierr = TaoLogConvergenceHistory(tao,asls->merit,ndpsi,0.0,tao->ksp_its);CHKERRQ(ierr);
158 ierr = TaoMonitor(tao,tao->niter,asls->merit,ndpsi,0.0,t);CHKERRQ(ierr);
159 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr);
160 if (TAO_CONTINUE_ITERATING != tao->reason) break;
161
162 /* Call general purpose update function */
163 if (tao->ops->update) {
164 ierr = (*tao->ops->update)(tao, tao->niter, tao->user_update);CHKERRQ(ierr);
165 }
166 tao->niter++;
167
168 /* We are going to solve a linear system of equations. We need to
169 set the tolerances for the solve so that we maintain an asymptotic
170 rate of convergence that is superlinear.
171 Note: these tolerances are for the reduced system. We really need
172 to make sure that the full system satisfies the full-space conditions.
173
174 This rule gives superlinear asymptotic convergence
175 asls->atol = min(0.5, asls->merit*sqrt(asls->merit));
176 asls->rtol = 0.0;
177
178 This rule gives quadratic asymptotic convergence
179 asls->atol = min(0.5, asls->merit*asls->merit);
180 asls->rtol = 0.0;
181
182 Calculate a free and fixed set of variables. The fixed set of
183 variables are those for the d_b is approximately equal to zero.
184 The definition of approximately changes as we approach the solution
185 to the problem.
186
187 No one rule is guaranteed to work in all cases. The following
188 definition is based on the norm of the Jacobian matrix. If the
189 norm is large, the tolerance becomes smaller. */
190 ierr = MatNorm(tao->jacobian,NORM_1,&asls->identifier);CHKERRQ(ierr);
191 asls->identifier = PetscMin(asls->merit, 1e-2) / (1 + asls->identifier);
192
193 ierr = VecSet(asls->t1,-asls->identifier);CHKERRQ(ierr);
194 ierr = VecSet(asls->t2, asls->identifier);CHKERRQ(ierr);
195
196 ierr = ISDestroy(&asls->fixed);CHKERRQ(ierr);
197 ierr = ISDestroy(&asls->free);CHKERRQ(ierr);
198 ierr = VecWhichBetweenOrEqual(asls->t1, asls->db, asls->t2, &asls->fixed);CHKERRQ(ierr);
199 ierr = ISComplementVec(asls->fixed,asls->t1, &asls->free);CHKERRQ(ierr);
200
201 ierr = ISGetSize(asls->fixed,&nf);CHKERRQ(ierr);
202 ierr = PetscInfo1(tao,"Number of fixed variables: %D\n", nf);CHKERRQ(ierr);
203
204 /* We now have our partition. Now calculate the direction in the
205 fixed variable space. */
206 ierr = TaoVecGetSubVec(asls->ff, asls->fixed, tao->subset_type, 0.0, &asls->r1);CHKERRQ(ierr);
207 ierr = TaoVecGetSubVec(asls->da, asls->fixed, tao->subset_type, 1.0, &asls->r2);CHKERRQ(ierr);
208 ierr = VecPointwiseDivide(asls->r1,asls->r1,asls->r2);CHKERRQ(ierr);
209 ierr = VecSet(tao->stepdirection,0.0);CHKERRQ(ierr);
210 ierr = VecISAXPY(tao->stepdirection, asls->fixed, 1.0,asls->r1);CHKERRQ(ierr);
211
212 /* Our direction in the Fixed Variable Set is fixed. Calculate the
213 information needed for the step in the Free Variable Set. To
214 do this, we need to know the diagonal perturbation and the
215 right hand side. */
216
217 ierr = TaoVecGetSubVec(asls->da, asls->free, tao->subset_type, 0.0, &asls->r1);CHKERRQ(ierr);
218 ierr = TaoVecGetSubVec(asls->ff, asls->free, tao->subset_type, 0.0, &asls->r2);CHKERRQ(ierr);
219 ierr = TaoVecGetSubVec(asls->db, asls->free, tao->subset_type, 1.0, &asls->r3);CHKERRQ(ierr);
220 ierr = VecPointwiseDivide(asls->r1,asls->r1, asls->r3);CHKERRQ(ierr);
221 ierr = VecPointwiseDivide(asls->r2,asls->r2, asls->r3);CHKERRQ(ierr);
222
223 /* r1 is the diagonal perturbation
224 r2 is the right hand side
225 r3 is no longer needed
226
227 Now need to modify r2 for our direction choice in the fixed
228 variable set: calculate t1 = J*d, take the reduced vector
229 of t1 and modify r2. */
230
231 ierr = MatMult(tao->jacobian, tao->stepdirection, asls->t1);CHKERRQ(ierr);
232 ierr = TaoVecGetSubVec(asls->t1,asls->free,tao->subset_type,0.0,&asls->r3);CHKERRQ(ierr);
233 ierr = VecAXPY(asls->r2, -1.0, asls->r3);CHKERRQ(ierr);
234
235 /* Calculate the reduced problem matrix and the direction */
236 ierr = TaoMatGetSubMat(tao->jacobian, asls->free, asls->w, tao->subset_type,&asls->J_sub);CHKERRQ(ierr);
237 if (tao->jacobian != tao->jacobian_pre) {
238 ierr = TaoMatGetSubMat(tao->jacobian_pre, asls->free, asls->w, tao->subset_type, &asls->Jpre_sub);CHKERRQ(ierr);
239 } else {
240 ierr = MatDestroy(&asls->Jpre_sub);CHKERRQ(ierr);
241 asls->Jpre_sub = asls->J_sub;
242 ierr = PetscObjectReference((PetscObject)(asls->Jpre_sub));CHKERRQ(ierr);
243 }
244 ierr = MatDiagonalSet(asls->J_sub, asls->r1,ADD_VALUES);CHKERRQ(ierr);
245 ierr = TaoVecGetSubVec(tao->stepdirection, asls->free, tao->subset_type, 0.0, &asls->dxfree);CHKERRQ(ierr);
246 ierr = VecSet(asls->dxfree, 0.0);CHKERRQ(ierr);
247
248 /* Calculate the reduced direction. (Really negative of Newton
249 direction. Therefore, rest of the code uses -d.) */
250 ierr = KSPReset(tao->ksp);CHKERRQ(ierr);
251 ierr = KSPSetOperators(tao->ksp, asls->J_sub, asls->Jpre_sub);CHKERRQ(ierr);
252 ierr = KSPSolve(tao->ksp, asls->r2, asls->dxfree);CHKERRQ(ierr);
253 ierr = KSPGetIterationNumber(tao->ksp,&tao->ksp_its);CHKERRQ(ierr);
254 tao->ksp_tot_its+=tao->ksp_its;
255
256 /* Add the direction in the free variables back into the real direction. */
257 ierr = VecISAXPY(tao->stepdirection, asls->free, 1.0,asls->dxfree);CHKERRQ(ierr);
258
259
260 /* Check the projected real direction for descent and if not, use the negative
261 gradient direction. */
262 ierr = VecCopy(tao->stepdirection, asls->w);CHKERRQ(ierr);
263 ierr = VecScale(asls->w, -1.0);CHKERRQ(ierr);
264 ierr = VecBoundGradientProjection(asls->w, tao->solution, tao->XL, tao->XU, asls->w);CHKERRQ(ierr);
265 ierr = VecNorm(asls->w, NORM_2, &normd);CHKERRQ(ierr);
266 ierr = VecDot(asls->w, asls->dpsi, &innerd);CHKERRQ(ierr);
267
268 if (innerd >= -asls->delta*PetscPowReal(normd, asls->rho)) {
269 ierr = PetscInfo1(tao,"Gradient direction: %5.4e.\n", (double)innerd);CHKERRQ(ierr);
270 ierr = PetscInfo1(tao, "Iteration %D: newton direction not descent\n", tao->niter);CHKERRQ(ierr);
271 ierr = VecCopy(asls->dpsi, tao->stepdirection);CHKERRQ(ierr);
272 ierr = VecDot(asls->dpsi, tao->stepdirection, &innerd);CHKERRQ(ierr);
273 }
274
275 ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
276 innerd = -innerd;
277
278 /* We now have a correct descent direction. Apply a linesearch to
279 find the new iterate. */
280 ierr = TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0);CHKERRQ(ierr);
281 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &psi,asls->dpsi, tao->stepdirection, &t, &ls_reason);CHKERRQ(ierr);
282 ierr = VecNorm(asls->dpsi, NORM_2, &ndpsi);CHKERRQ(ierr);
283 }
284 PetscFunctionReturn(0);
285 }
286
287 /* ---------------------------------------------------------- */
288 /*MC
289 TAOASFLS - Active-set feasible linesearch algorithm for solving
290 complementarity constraints
291
292 Options Database Keys:
293 + -tao_ssls_delta - descent test fraction
294 - -tao_ssls_rho - descent test power
295
296 Level: beginner
297 M*/
TaoCreate_ASFLS(Tao tao)298 PETSC_EXTERN PetscErrorCode TaoCreate_ASFLS(Tao tao)
299 {
300 TAO_SSLS *asls;
301 PetscErrorCode ierr;
302 const char *armijo_type = TAOLINESEARCHARMIJO;
303
304 PetscFunctionBegin;
305 ierr = PetscNewLog(tao,&asls);CHKERRQ(ierr);
306 tao->data = (void*)asls;
307 tao->ops->solve = TaoSolve_ASFLS;
308 tao->ops->setup = TaoSetUp_ASFLS;
309 tao->ops->view = TaoView_SSLS;
310 tao->ops->setfromoptions = TaoSetFromOptions_SSLS;
311 tao->ops->destroy = TaoDestroy_ASFLS;
312 tao->subset_type = TAO_SUBSET_SUBVEC;
313 asls->delta = 1e-10;
314 asls->rho = 2.1;
315 asls->fixed = NULL;
316 asls->free = NULL;
317 asls->J_sub = NULL;
318 asls->Jpre_sub = NULL;
319 asls->w = NULL;
320 asls->r1 = NULL;
321 asls->r2 = NULL;
322 asls->r3 = NULL;
323 asls->t1 = NULL;
324 asls->t2 = NULL;
325 asls->dxfree = NULL;
326 asls->identifier = 1e-5;
327
328 ierr = TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);CHKERRQ(ierr);
329 ierr = PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);CHKERRQ(ierr);
330 ierr = TaoLineSearchSetType(tao->linesearch, armijo_type);CHKERRQ(ierr);
331 ierr = TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);CHKERRQ(ierr);
332 ierr = TaoLineSearchSetFromOptions(tao->linesearch);CHKERRQ(ierr);
333
334 ierr = KSPCreate(((PetscObject)tao)->comm, &tao->ksp);CHKERRQ(ierr);
335 ierr = PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);CHKERRQ(ierr);
336 ierr = KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix);CHKERRQ(ierr);
337 ierr = KSPSetFromOptions(tao->ksp);CHKERRQ(ierr);
338
339 /* Override default settings (unless already changed) */
340 if (!tao->max_it_changed) tao->max_it = 2000;
341 if (!tao->max_funcs_changed) tao->max_funcs = 4000;
342 if (!tao->gttol_changed) tao->gttol = 0;
343 if (!tao->grtol_changed) tao->grtol = 0;
344 #if defined(PETSC_USE_REAL_SINGLE)
345 if (!tao->gatol_changed) tao->gatol = 1.0e-6;
346 if (!tao->fmin_changed) tao->fmin = 1.0e-4;
347 #else
348 if (!tao->gatol_changed) tao->gatol = 1.0e-16;
349 if (!tao->fmin_changed) tao->fmin = 1.0e-8;
350 #endif
351 PetscFunctionReturn(0);
352 }
353