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