Actual source code: asfls.c
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
14: Many other reformulations are possible and combinations of
15: feasible/infeasible and linesearch/trust region are possible.
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.
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.
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.
36: The theory for the feasible version follows from the feasible descent
37: algorithm framework.
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: */
55: static PetscErrorCode TaoSetUp_ASFLS(Tao tao)
56: {
57: TAO_SSLS *asls = (TAO_SSLS *)tao->data;
59: PetscFunctionBegin;
60: PetscCall(VecDuplicate(tao->solution, &tao->gradient));
61: PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
62: PetscCall(VecDuplicate(tao->solution, &asls->ff));
63: PetscCall(VecDuplicate(tao->solution, &asls->dpsi));
64: PetscCall(VecDuplicate(tao->solution, &asls->da));
65: PetscCall(VecDuplicate(tao->solution, &asls->db));
66: PetscCall(VecDuplicate(tao->solution, &asls->t1));
67: PetscCall(VecDuplicate(tao->solution, &asls->t2));
68: PetscCall(VecDuplicate(tao->solution, &asls->w));
69: asls->fixed = NULL;
70: asls->free = NULL;
71: asls->J_sub = NULL;
72: asls->Jpre_sub = NULL;
73: asls->r1 = NULL;
74: asls->r2 = NULL;
75: asls->r3 = NULL;
76: asls->dxfree = NULL;
77: PetscFunctionReturn(PETSC_SUCCESS);
78: }
80: static PetscErrorCode Tao_ASLS_FunctionGradient(TaoLineSearch ls, Vec X, PetscReal *fcn, Vec G, void *ptr)
81: {
82: Tao tao = (Tao)ptr;
83: TAO_SSLS *asls = (TAO_SSLS *)tao->data;
85: PetscFunctionBegin;
86: PetscCall(TaoComputeConstraints(tao, X, tao->constraints));
87: PetscCall(VecFischer(X, tao->constraints, tao->XL, tao->XU, asls->ff));
88: PetscCall(VecNorm(asls->ff, NORM_2, &asls->merit));
89: *fcn = 0.5 * asls->merit * asls->merit;
90: PetscCall(TaoComputeJacobian(tao, tao->solution, tao->jacobian, tao->jacobian_pre));
92: PetscCall(MatDFischer(tao->jacobian, tao->solution, tao->constraints, tao->XL, tao->XU, asls->t1, asls->t2, asls->da, asls->db));
93: PetscCall(VecPointwiseMult(asls->t1, asls->ff, asls->db));
94: PetscCall(MatMultTranspose(tao->jacobian, asls->t1, G));
95: PetscCall(VecPointwiseMult(asls->t1, asls->ff, asls->da));
96: PetscCall(VecAXPY(G, 1.0, asls->t1));
97: PetscFunctionReturn(PETSC_SUCCESS);
98: }
100: static PetscErrorCode TaoDestroy_ASFLS(Tao tao)
101: {
102: TAO_SSLS *ssls = (TAO_SSLS *)tao->data;
104: PetscFunctionBegin;
105: PetscCall(VecDestroy(&ssls->ff));
106: PetscCall(VecDestroy(&ssls->dpsi));
107: PetscCall(VecDestroy(&ssls->da));
108: PetscCall(VecDestroy(&ssls->db));
109: PetscCall(VecDestroy(&ssls->w));
110: PetscCall(VecDestroy(&ssls->t1));
111: PetscCall(VecDestroy(&ssls->t2));
112: PetscCall(VecDestroy(&ssls->r1));
113: PetscCall(VecDestroy(&ssls->r2));
114: PetscCall(VecDestroy(&ssls->r3));
115: PetscCall(VecDestroy(&ssls->dxfree));
116: PetscCall(MatDestroy(&ssls->J_sub));
117: PetscCall(MatDestroy(&ssls->Jpre_sub));
118: PetscCall(ISDestroy(&ssls->fixed));
119: PetscCall(ISDestroy(&ssls->free));
120: PetscCall(KSPDestroy(&tao->ksp));
121: PetscCall(PetscFree(tao->data));
122: PetscFunctionReturn(PETSC_SUCCESS);
123: }
125: static PetscErrorCode TaoSolve_ASFLS(Tao tao)
126: {
127: TAO_SSLS *asls = (TAO_SSLS *)tao->data;
128: PetscReal psi, ndpsi, normd, innerd, t = 0;
129: PetscInt nf;
130: TaoLineSearchConvergedReason ls_reason;
132: PetscFunctionBegin;
133: /* Assume that Setup has been called!
134: Set the structure for the Jacobian and create a linear solver. */
136: PetscCall(TaoComputeVariableBounds(tao));
137: PetscCall(TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch, Tao_ASLS_FunctionGradient, tao));
138: PetscCall(TaoLineSearchSetObjectiveRoutine(tao->linesearch, Tao_SSLS_Function, tao));
139: PetscCall(TaoLineSearchSetVariableBounds(tao->linesearch, tao->XL, tao->XU));
141: PetscCall(VecMedian(tao->XL, tao->solution, tao->XU, tao->solution));
143: /* Calculate the function value and fischer function value at the
144: current iterate */
145: PetscCall(TaoLineSearchComputeObjectiveAndGradient(tao->linesearch, tao->solution, &psi, asls->dpsi));
146: PetscCall(VecNorm(asls->dpsi, NORM_2, &ndpsi));
148: tao->reason = TAO_CONTINUE_ITERATING;
149: while (1) {
150: /* Check the converged criteria */
151: PetscCall(PetscInfo(tao, "iter %" PetscInt_FMT ", merit: %g, ||dpsi||: %g\n", tao->niter, (double)asls->merit, (double)ndpsi));
152: PetscCall(TaoLogConvergenceHistory(tao, asls->merit, ndpsi, 0.0, tao->ksp_its));
153: PetscCall(TaoMonitor(tao, tao->niter, asls->merit, ndpsi, 0.0, t));
154: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
155: if (TAO_CONTINUE_ITERATING != tao->reason) break;
157: /* Call general purpose update function */
158: PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
159: tao->niter++;
161: /* We are going to solve a linear system of equations. We need to
162: set the tolerances for the solve so that we maintain an asymptotic
163: rate of convergence that is superlinear.
164: Note: these tolerances are for the reduced system. We really need
165: to make sure that the full system satisfies the full-space conditions.
167: This rule gives superlinear asymptotic convergence
168: asls->atol = min(0.5, asls->merit*sqrt(asls->merit));
169: asls->rtol = 0.0;
171: This rule gives quadratic asymptotic convergence
172: asls->atol = min(0.5, asls->merit*asls->merit);
173: asls->rtol = 0.0;
175: Calculate a free and fixed set of variables. The fixed set of
176: variables are those for the d_b is approximately equal to zero.
177: The definition of approximately changes as we approach the solution
178: to the problem.
180: No one rule is guaranteed to work in all cases. The following
181: definition is based on the norm of the Jacobian matrix. If the
182: norm is large, the tolerance becomes smaller. */
183: PetscCall(MatNorm(tao->jacobian, NORM_1, &asls->identifier));
184: asls->identifier = PetscMin(asls->merit, 1e-2) / (1 + asls->identifier);
186: PetscCall(VecSet(asls->t1, -asls->identifier));
187: PetscCall(VecSet(asls->t2, asls->identifier));
189: PetscCall(ISDestroy(&asls->fixed));
190: PetscCall(ISDestroy(&asls->free));
191: PetscCall(VecWhichBetweenOrEqual(asls->t1, asls->db, asls->t2, &asls->fixed));
192: PetscCall(ISComplementVec(asls->fixed, asls->t1, &asls->free));
194: PetscCall(ISGetSize(asls->fixed, &nf));
195: PetscCall(PetscInfo(tao, "Number of fixed variables: %" PetscInt_FMT "\n", nf));
197: /* We now have our partition. Now calculate the direction in the
198: fixed variable space. */
199: PetscCall(TaoVecGetSubVec(asls->ff, asls->fixed, tao->subset_type, 0.0, &asls->r1));
200: PetscCall(TaoVecGetSubVec(asls->da, asls->fixed, tao->subset_type, 1.0, &asls->r2));
201: PetscCall(VecPointwiseDivide(asls->r1, asls->r1, asls->r2));
202: PetscCall(VecSet(tao->stepdirection, 0.0));
203: PetscCall(VecISAXPY(tao->stepdirection, asls->fixed, 1.0, asls->r1));
205: /* Our direction in the Fixed Variable Set is fixed. Calculate the
206: information needed for the step in the Free Variable Set. To
207: do this, we need to know the diagonal perturbation and the
208: right hand side. */
210: PetscCall(TaoVecGetSubVec(asls->da, asls->free, tao->subset_type, 0.0, &asls->r1));
211: PetscCall(TaoVecGetSubVec(asls->ff, asls->free, tao->subset_type, 0.0, &asls->r2));
212: PetscCall(TaoVecGetSubVec(asls->db, asls->free, tao->subset_type, 1.0, &asls->r3));
213: PetscCall(VecPointwiseDivide(asls->r1, asls->r1, asls->r3));
214: PetscCall(VecPointwiseDivide(asls->r2, asls->r2, asls->r3));
216: /* r1 is the diagonal perturbation
217: r2 is the right hand side
218: r3 is no longer needed
220: Now need to modify r2 for our direction choice in the fixed
221: variable set: calculate t1 = J*d, take the reduced vector
222: of t1 and modify r2. */
224: PetscCall(MatMult(tao->jacobian, tao->stepdirection, asls->t1));
225: PetscCall(TaoVecGetSubVec(asls->t1, asls->free, tao->subset_type, 0.0, &asls->r3));
226: PetscCall(VecAXPY(asls->r2, -1.0, asls->r3));
228: /* Calculate the reduced problem matrix and the direction */
229: PetscCall(TaoMatGetSubMat(tao->jacobian, asls->free, asls->w, tao->subset_type, &asls->J_sub));
230: if (tao->jacobian != tao->jacobian_pre) {
231: PetscCall(TaoMatGetSubMat(tao->jacobian_pre, asls->free, asls->w, tao->subset_type, &asls->Jpre_sub));
232: } else {
233: PetscCall(MatDestroy(&asls->Jpre_sub));
234: asls->Jpre_sub = asls->J_sub;
235: PetscCall(PetscObjectReference((PetscObject)(asls->Jpre_sub)));
236: }
237: PetscCall(MatDiagonalSet(asls->J_sub, asls->r1, ADD_VALUES));
238: PetscCall(TaoVecGetSubVec(tao->stepdirection, asls->free, tao->subset_type, 0.0, &asls->dxfree));
239: PetscCall(VecSet(asls->dxfree, 0.0));
241: /* Calculate the reduced direction. (Really negative of Newton
242: direction. Therefore, rest of the code uses -d.) */
243: PetscCall(KSPReset(tao->ksp));
244: PetscCall(KSPSetOperators(tao->ksp, asls->J_sub, asls->Jpre_sub));
245: PetscCall(KSPSolve(tao->ksp, asls->r2, asls->dxfree));
246: PetscCall(KSPGetIterationNumber(tao->ksp, &tao->ksp_its));
247: tao->ksp_tot_its += tao->ksp_its;
249: /* Add the direction in the free variables back into the real direction. */
250: PetscCall(VecISAXPY(tao->stepdirection, asls->free, 1.0, asls->dxfree));
252: /* Check the projected real direction for descent and if not, use the negative
253: gradient direction. */
254: PetscCall(VecCopy(tao->stepdirection, asls->w));
255: PetscCall(VecScale(asls->w, -1.0));
256: PetscCall(VecBoundGradientProjection(asls->w, tao->solution, tao->XL, tao->XU, asls->w));
257: PetscCall(VecNorm(asls->w, NORM_2, &normd));
258: PetscCall(VecDot(asls->w, asls->dpsi, &innerd));
260: if (innerd >= -asls->delta * PetscPowReal(normd, asls->rho)) {
261: PetscCall(PetscInfo(tao, "Gradient direction: %5.4e.\n", (double)innerd));
262: PetscCall(PetscInfo(tao, "Iteration %" PetscInt_FMT ": newton direction not descent\n", tao->niter));
263: PetscCall(VecCopy(asls->dpsi, tao->stepdirection));
264: PetscCall(VecDot(asls->dpsi, tao->stepdirection, &innerd));
265: }
267: PetscCall(VecScale(tao->stepdirection, -1.0));
268: innerd = -innerd;
270: /* We now have a correct descent direction. Apply a linesearch to
271: find the new iterate. */
272: PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
273: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &psi, asls->dpsi, tao->stepdirection, &t, &ls_reason));
274: PetscCall(VecNorm(asls->dpsi, NORM_2, &ndpsi));
275: }
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }
279: /* ---------------------------------------------------------- */
280: /*MC
281: TAOASFLS - Active-set feasible linesearch algorithm for solving
282: complementarity constraints
284: Options Database Keys:
285: + -tao_ssls_delta - descent test fraction
286: - -tao_ssls_rho - descent test power
288: Level: beginner
289: M*/
290: PETSC_EXTERN PetscErrorCode TaoCreate_ASFLS(Tao tao)
291: {
292: TAO_SSLS *asls;
293: const char *armijo_type = TAOLINESEARCHARMIJO;
295: PetscFunctionBegin;
296: PetscCall(PetscNew(&asls));
297: tao->data = (void *)asls;
298: tao->ops->solve = TaoSolve_ASFLS;
299: tao->ops->setup = TaoSetUp_ASFLS;
300: tao->ops->view = TaoView_SSLS;
301: tao->ops->setfromoptions = TaoSetFromOptions_SSLS;
302: tao->ops->destroy = TaoDestroy_ASFLS;
303: tao->subset_type = TAO_SUBSET_SUBVEC;
304: asls->delta = 1e-10;
305: asls->rho = 2.1;
306: asls->fixed = NULL;
307: asls->free = NULL;
308: asls->J_sub = NULL;
309: asls->Jpre_sub = NULL;
310: asls->w = NULL;
311: asls->r1 = NULL;
312: asls->r2 = NULL;
313: asls->r3 = NULL;
314: asls->t1 = NULL;
315: asls->t2 = NULL;
316: asls->dxfree = NULL;
317: asls->identifier = 1e-5;
319: PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
320: PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
321: PetscCall(TaoLineSearchSetType(tao->linesearch, armijo_type));
322: PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix));
323: PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
325: PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp));
326: PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1));
327: PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix));
328: PetscCall(KSPSetFromOptions(tao->ksp));
330: /* Override default settings (unless already changed) */
331: if (!tao->max_it_changed) tao->max_it = 2000;
332: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
333: if (!tao->gttol_changed) tao->gttol = 0;
334: if (!tao->grtol_changed) tao->grtol = 0;
335: #if defined(PETSC_USE_REAL_SINGLE)
336: if (!tao->gatol_changed) tao->gatol = 1.0e-6;
337: if (!tao->fmin_changed) tao->fmin = 1.0e-4;
338: #else
339: if (!tao->gatol_changed) tao->gatol = 1.0e-16;
340: if (!tao->fmin_changed) tao->fmin = 1.0e-8;
341: #endif
342: PetscFunctionReturn(PETSC_SUCCESS);
343: }