Actual source code: asfls.c
1: #include "src/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 407-439, 1996.
45: Ferris, Kanzow, Munson, "Feasible Descent Algorithms for Mixed
46: Complementarity Problems," Mathematical Programming, 86,
47: pages 475-497, 1999.
48: Fischer, "A Special Newton-type Optimization Method," Optimization,
49: 24, pages 269-284, 1992
50: Munson, Facchinei, Ferris, Fischer, Kanzow, "The Semismooth Algorithm
51: for Large Scale Complementarity Problems," Technical Report 99-06,
52: University of Wisconsin - Madison, 1999.
53: */
58: PetscErrorCode TaoSetUp_ASFLS(TaoSolver tao)
59: {
60: TAO_SSLS *asls = (TAO_SSLS *)tao->data;
65: VecDuplicate(tao->solution,&tao->gradient);
66: VecDuplicate(tao->solution,&tao->stepdirection);
67: VecDuplicate(tao->solution,&asls->ff);
68: VecDuplicate(tao->solution,&asls->dpsi);
69: VecDuplicate(tao->solution,&asls->da);
70: VecDuplicate(tao->solution,&asls->db);
71: VecDuplicate(tao->solution,&asls->t1);
72: VecDuplicate(tao->solution,&asls->t2);
73: VecDuplicate(tao->solution, &asls->w);
74: asls->fixed = PETSC_NULL;
75: asls->free = PETSC_NULL;
76: asls->J_sub = PETSC_NULL;
77: asls->Jpre_sub = PETSC_NULL;
78: asls->r1 = PETSC_NULL;
79: asls->r2 = PETSC_NULL;
80: asls->r3 = PETSC_NULL;
81: asls->dxfree = PETSC_NULL;
83: return(0);
84: }
88: static PetscErrorCode Tao_ASLS_FunctionGradient(TaoLineSearch ls, Vec X, PetscReal *fcn, Vec G, void *ptr)
89: {
90: TaoSolver tao = (TaoSolver)ptr;
91: TAO_SSLS *asls = (TAO_SSLS *)tao->data;
96: TaoComputeConstraints(tao, X, tao->constraints);
97: VecFischer(X,tao->constraints,tao->XL,tao->XU,asls->ff);
98: VecNorm(asls->ff,NORM_2,&asls->merit);
99: *fcn = 0.5*asls->merit*asls->merit;
101: TaoComputeJacobian(tao, tao->solution, &tao->jacobian, &tao->jacobian_pre, &asls->matflag);
102:
103: D_Fischer(tao->jacobian, tao->solution, tao->constraints,
104: tao->XL, tao->XU, asls->t1, asls->t2,
105: asls->da, asls->db);
106: VecPointwiseMult(asls->t1, asls->ff, asls->db);
107: MatMultTranspose(tao->jacobian,asls->t1,G);
108: VecPointwiseMult(asls->t1, asls->ff, asls->da);
109: VecAXPY(G,1.0,asls->t1);
110: return(0);
111: }
115: static PetscErrorCode TaoDestroy_ASFLS(TaoSolver tao)
116: {
117: TAO_SSLS *ssls = (TAO_SSLS *)tao->data;
122: VecDestroy(&ssls->ff);
123: VecDestroy(&ssls->dpsi);
124: VecDestroy(&ssls->da);
125: VecDestroy(&ssls->db);
126: VecDestroy(&ssls->w);
127: VecDestroy(&ssls->t1);
128: VecDestroy(&ssls->t2);
129: VecDestroy(&ssls->r1);
130: VecDestroy(&ssls->r2);
131: VecDestroy(&ssls->r3);
132: VecDestroy(&ssls->dxfree);
133: MatDestroy(&ssls->J_sub);
134: MatDestroy(&ssls->Jpre_sub);
135: ISDestroy(&ssls->fixed);
136: ISDestroy(&ssls->free);
137: PetscFree(tao->data);
138:
139: tao->data = PETSC_NULL;
140: return(0);
141:
142: }
145: static PetscErrorCode TaoSolve_ASFLS(TaoSolver tao)
146: {
147: TAO_SSLS *asls = (TAO_SSLS *)tao->data;
148: PetscReal psi,ndpsi, normd, innerd, t=0;
149: PetscInt iter=0, nf;
151: TaoSolverTerminationReason reason;
152: TaoLineSearchTerminationReason ls_reason;
156: /* Assume that Setup has been called!
157: Set the structure for the Jacobian and create a linear solver. */
159: TaoComputeVariableBounds(tao);
160: TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,Tao_ASLS_FunctionGradient,tao);
161: TaoLineSearchSetObjectiveRoutine(tao->linesearch,Tao_SSLS_Function,tao);
162: TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);
164: VecMedian(tao->XL, tao->solution, tao->XU, tao->solution);
166: /* Calculate the function value and fischer function value at the
167: current iterate */
168: TaoLineSearchComputeObjectiveAndGradient(tao->linesearch,tao->solution,&psi,asls->dpsi);
169: VecNorm(asls->dpsi,NORM_2,&ndpsi);
171: while (1) {
173: /* Check the termination criteria */
174: PetscInfo3(tao,"iter %D, merit: %G, ||dpsi||: %G\n",iter, asls->merit, ndpsi);
175: TaoMonitor(tao, iter++, asls->merit, ndpsi, 0.0, t, &reason);
176: if (TAO_CONTINUE_ITERATING != reason) break;
178: /* We are going to solve a linear system of equations. We need to
179: set the tolerances for the solve so that we maintain an asymptotic
180: rate of convergence that is superlinear.
181: Note: these tolerances are for the reduced system. We really need
182: to make sure that the full system satisfies the full-space conditions.
183:
184: This rule gives superlinear asymptotic convergence
185: asls->atol = min(0.5, asls->merit*sqrt(asls->merit));
186: asls->rtol = 0.0;
188: This rule gives quadratic asymptotic convergence
189: asls->atol = min(0.5, asls->merit*asls->merit);
190: asls->rtol = 0.0;
192: Calculate a free and fixed set of variables. The fixed set of
193: variables are those for the d_b is approximately equal to zero.
194: The definition of approximately changes as we approach the solution
195: to the problem.
197: No one rule is guaranteed to work in all cases. The following
198: definition is based on the norm of the Jacobian matrix. If the
199: norm is large, the tolerance becomes smaller. */
200: MatNorm(tao->jacobian,NORM_1,&asls->identifier);
201: asls->identifier = PetscMin(asls->merit, 1e-2) / (1 + asls->identifier);
203: VecSet(asls->t1,-asls->identifier);
204: VecSet(asls->t2, asls->identifier);
206: ISDestroy(&asls->fixed);
207: ISDestroy(&asls->free);
208: VecWhichBetweenOrEqual(asls->t1, asls->db, asls->t2, &asls->fixed);
209: ISCreateComplement(asls->fixed,asls->t1, &asls->free);
211: ISGetSize(asls->fixed,&nf);
212: PetscInfo1(tao,"Number of fixed variables: %d\n", nf);
214: /* We now have our partition. Now calculate the direction in the
215: fixed variable space. */
216: VecGetSubVec(asls->ff, asls->fixed, tao->subset_type, 0.0, &asls->r1);
217: VecGetSubVec(asls->da, asls->fixed, tao->subset_type, 1.0, &asls->r2);
218: VecPointwiseDivide(asls->r1,asls->r1,asls->r2);
219: VecSet(tao->stepdirection,0.0);
220: VecReducedXPY(tao->stepdirection,asls->r1, asls->fixed);
223: /* Our direction in the Fixed Variable Set is fixed. Calculate the
224: information needed for the step in the Free Variable Set. To
225: do this, we need to know the diagonal perturbation and the
226: right hand side. */
228: VecGetSubVec(asls->da, asls->free, tao->subset_type, 0.0, &asls->r1);
229: VecGetSubVec(asls->ff, asls->free, tao->subset_type, 0.0, &asls->r2);
230: VecGetSubVec(asls->db, asls->free, tao->subset_type, 1.0, &asls->r3);
231: VecPointwiseDivide(asls->r1,asls->r1, asls->r3);
232: VecPointwiseDivide(asls->r2,asls->r2, asls->r3);
234: /* r1 is the diagonal perturbation
235: r2 is the right hand side
236: r3 is no longer needed
238: Now need to modify r2 for our direction choice in the fixed
239: variable set: calculate t1 = J*d, take the reduced vector
240: of t1 and modify r2. */
242: MatMult(tao->jacobian, tao->stepdirection, asls->t1);
243: VecGetSubVec(asls->t1,asls->free,tao->subset_type,0.0,&asls->r3);
244: VecAXPY(asls->r2, -1.0, asls->r3);
246: /* Calculate the reduced problem matrix and the direction */
247: MatGetSubMat(tao->jacobian, asls->free, asls->w, tao->subset_type,&asls->J_sub);
248: if (tao->jacobian != tao->jacobian_pre) {
249: MatGetSubMat(tao->jacobian_pre, asls->free, asls->w, tao->subset_type, &asls->Jpre_sub);
250: } else {
251: MatDestroy(&asls->Jpre_sub);
252: asls->Jpre_sub = asls->J_sub;
253: PetscObjectReference((PetscObject)(asls->Jpre_sub));
254: }
255: MatDiagonalSet(asls->J_sub, asls->r1,ADD_VALUES);
256: VecGetSubVec(tao->stepdirection, asls->free, tao->subset_type, 0.0, &asls->dxfree);
257: VecSet(asls->dxfree, 0.0);
259: /* Calculate the reduced direction. (Really negative of Newton
260: direction. Therefore, rest of the code uses -d.) */
261: KSPReset(tao->ksp);
262: KSPSetOperators(tao->ksp, asls->J_sub, asls->Jpre_sub, asls->matflag);
263: KSPSolve(tao->ksp, asls->r2, asls->dxfree);
265: /* Add the direction in the free variables back into the real direction. */
266: VecReducedXPY(tao->stepdirection, asls->dxfree, asls->free);
269: /* Check the projected real direction for descent and if not, use the negative
270: gradient direction. */
271: VecCopy(tao->stepdirection, asls->w);
272: VecScale(asls->w, -1.0);
273: VecBoundGradientProjection(asls->w, tao->solution, tao->XL, tao->XU, asls->w);
274: VecNorm(asls->w, NORM_2, &normd);
275: VecDot(asls->w, asls->dpsi, &innerd);
277: if (innerd >= -asls->delta*pow(normd, asls->rho)) {
278: PetscInfo1(tao,"Gradient direction: %5.4e.\n", innerd);
279: PetscInfo1(tao, "Iteration %d: newton direction not descent\n", iter);
280: VecCopy(asls->dpsi, tao->stepdirection);
281: VecDot(asls->dpsi, tao->stepdirection, &innerd);
282: }
284: VecScale(tao->stepdirection, -1.0);
285: innerd = -innerd;
287: /* We now have a correct descent direction. Apply a linesearch to
288: find the new iterate. */
289: TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0);
290: TaoLineSearchApply(tao->linesearch, tao->solution, &psi,
291: asls->dpsi, tao->stepdirection, &t, &ls_reason);
292: VecNorm(asls->dpsi, NORM_2, &ndpsi);
293: }
295: return(0);
296: }
298: /* ---------------------------------------------------------- */
302: PetscErrorCode TaoCreate_ASFLS(TaoSolver tao)
303: {
304: TAO_SSLS *asls;
305: PetscErrorCode ierr;
306: const char *armijo_type = TAOLINESEARCH_ARMIJO;
309: PetscNewLog(tao,TAO_SSLS,&asls);
310: tao->data = (void*)asls;
311: tao->ops->solve = TaoSolve_ASFLS;
312: tao->ops->setup = TaoSetUp_ASFLS;
313: tao->ops->view = TaoView_SSLS;
314: tao->ops->setfromoptions = TaoSetFromOptions_SSLS;
315: tao->ops->destroy = TaoDestroy_ASFLS;
316: tao->subset_type = TAO_SUBSET_SUBVEC;
317: asls->delta = 1e-10;
318: asls->rho = 2.1;
319: asls->fixed = PETSC_NULL;
320: asls->free = PETSC_NULL;
321: asls->J_sub = PETSC_NULL;
322: asls->Jpre_sub = PETSC_NULL;
323: asls->w = PETSC_NULL;
324: asls->r1 = PETSC_NULL;
325: asls->r2 = PETSC_NULL;
326: asls->r3 = PETSC_NULL;
327: asls->t1 = PETSC_NULL;
328: asls->t2 = PETSC_NULL;
329: asls->dxfree = PETSC_NULL;
330:
331: asls->identifier = 1e-5;
333: TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
334: TaoLineSearchSetType(tao->linesearch, armijo_type);
335: TaoLineSearchSetFromOptions(tao->linesearch);
337: KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
338: KSPSetFromOptions(tao->ksp);
339: tao->max_it = 2000;
340: tao->max_funcs = 4000;
341: tao->fatol = 0;
342: tao->frtol = 0;
343: tao->gttol = 0;
344: tao->grtol = 0;
345: tao->gatol = 1.0e-16;
346: tao->fmin = 1.0e-8;
348: return(0);
349: }