Actual source code: asils.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_ASILS(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: asls->fixed = PETSC_NULL;
74: asls->free = PETSC_NULL;
75: asls->J_sub = PETSC_NULL;
76: asls->Jpre_sub = PETSC_NULL;
77: asls->w = 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_ASILS(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_ASILS(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);
164: /* Calculate the function value and fischer function value at the
165: current iterate */
166: TaoLineSearchComputeObjectiveAndGradient(tao->linesearch,tao->solution,&psi,asls->dpsi);
167: VecNorm(asls->dpsi,NORM_2,&ndpsi);
169: while (1) {
171: /* Check the termination criteria */
172: PetscInfo3(tao,"iter %D, merit: %G, ||dpsi||: %G\n",iter, asls->merit, ndpsi);
173: TaoMonitor(tao, iter++, asls->merit, ndpsi, 0.0, t, &reason);
174: if (TAO_CONTINUE_ITERATING != reason) break;
176: /* We are going to solve a linear system of equations. We need to
177: set the tolerances for the solve so that we maintain an asymptotic
178: rate of convergence that is superlinear.
179: Note: these tolerances are for the reduced system. We really need
180: to make sure that the full system satisfies the full-space conditions.
181:
182: This rule gives superlinear asymptotic convergence
183: asls->atol = min(0.5, asls->merit*sqrt(asls->merit));
184: asls->rtol = 0.0;
186: This rule gives quadratic asymptotic convergence
187: asls->atol = min(0.5, asls->merit*asls->merit);
188: asls->rtol = 0.0;
190: Calculate a free and fixed set of variables. The fixed set of
191: variables are those for the d_b is approximately equal to zero.
192: The definition of approximately changes as we approach the solution
193: to the problem.
195: No one rule is guaranteed to work in all cases. The following
196: definition is based on the norm of the Jacobian matrix. If the
197: norm is large, the tolerance becomes smaller. */
198: MatNorm(tao->jacobian,NORM_1,&asls->identifier);
199: asls->identifier = PetscMin(asls->merit, 1e-2) / (1 + asls->identifier);
201: VecSet(asls->t1,-asls->identifier);
202: VecSet(asls->t2, asls->identifier);
204: ISDestroy(&asls->fixed);
205: ISDestroy(&asls->free);
206: VecWhichBetweenOrEqual(asls->t1, asls->db, asls->t2, &asls->fixed);
207: ISCreateComplement(asls->fixed,asls->t1, &asls->free);
209: ISGetSize(asls->fixed,&nf);
210: PetscInfo1(tao,"Number of fixed variables: %d\n", nf);
212: /* We now have our partition. Now calculate the direction in the
213: fixed variable space. */
214: VecGetSubVec(asls->ff, asls->fixed, tao->subset_type, 0.0, &asls->r1);
215: VecGetSubVec(asls->da, asls->fixed, tao->subset_type, 1.0, &asls->r2);
216: VecPointwiseDivide(asls->r1,asls->r1,asls->r2);
217: VecSet(tao->stepdirection,0.0);
218: VecReducedXPY(tao->stepdirection,asls->r1, asls->fixed);
221: /* Our direction in the Fixed Variable Set is fixed. Calculate the
222: information needed for the step in the Free Variable Set. To
223: do this, we need to know the diagonal perturbation and the
224: right hand side. */
226: VecGetSubVec(asls->da, asls->free, tao->subset_type, 0.0, &asls->r1);
227: VecGetSubVec(asls->ff, asls->free, tao->subset_type, 0.0, &asls->r2);
228: VecGetSubVec(asls->db, asls->free, tao->subset_type, 1.0, &asls->r3);
229: VecPointwiseDivide(asls->r1,asls->r1, asls->r3);
230: VecPointwiseDivide(asls->r2,asls->r2, asls->r3);
232: /* r1 is the diagonal perturbation
233: r2 is the right hand side
234: r3 is no longer needed
236: Now need to modify r2 for our direction choice in the fixed
237: variable set: calculate t1 = J*d, take the reduced vector
238: of t1 and modify r2. */
240: MatMult(tao->jacobian, tao->stepdirection, asls->t1);
241: VecGetSubVec(asls->t1,asls->free,tao->subset_type,0.0,&asls->r3);
242: VecAXPY(asls->r2, -1.0, asls->r3);
244: /* Calculate the reduced problem matrix and the direction */
245: if (!asls->w && (tao->subset_type == TAO_SUBSET_MASK
246: || tao->subset_type == TAO_SUBSET_MATRIXFREE)) {
247: VecDuplicate(tao->solution, &asls->w);
248: }
249: MatGetSubMat(tao->jacobian, asls->free, asls->w, tao->subset_type,&asls->J_sub);
250: if (tao->jacobian != tao->jacobian_pre) {
251: MatGetSubMat(tao->jacobian_pre, asls->free, asls->w, tao->subset_type, &asls->Jpre_sub);
252: } else {
253: MatDestroy(&asls->Jpre_sub);
254: asls->Jpre_sub = asls->J_sub;
255: PetscObjectReference((PetscObject)(asls->Jpre_sub));
256: }
257: MatDiagonalSet(asls->J_sub, asls->r1,ADD_VALUES);
258: VecGetSubVec(tao->stepdirection, asls->free, tao->subset_type, 0.0, &asls->dxfree);
259: VecSet(asls->dxfree, 0.0);
261: /* Calculate the reduced direction. (Really negative of Newton
262: direction. Therefore, rest of the code uses -d.) */
263: KSPReset(tao->ksp);
264: KSPSetOperators(tao->ksp, asls->J_sub, asls->Jpre_sub, asls->matflag);
265: KSPSolve(tao->ksp, asls->r2, asls->dxfree);
267: /* Add the direction in the free variables back into the real direction. */
268: VecReducedXPY(tao->stepdirection, asls->dxfree, asls->free);
271: /* Check the real direction for descent and if not, use the negative
272: gradient direction. */
273: VecNorm(tao->stepdirection, NORM_2, &normd);
274: VecDot(tao->stepdirection, asls->dpsi, &innerd);
276: if (innerd <= asls->delta*pow(normd, asls->rho)) {
277: PetscInfo1(tao,"Gradient direction: %5.4e.\n", innerd);
278: PetscInfo1(tao, "Iteration %d: newton direction not descent\n", iter);
279: VecCopy(asls->dpsi, tao->stepdirection);
280: VecDot(asls->dpsi, tao->stepdirection, &innerd);
281: }
283: VecScale(tao->stepdirection, -1.0);
284: innerd = -innerd;
286: /* We now have a correct descent direction. Apply a linesearch to
287: find the new iterate. */
288: TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0);
289: TaoLineSearchApply(tao->linesearch, tao->solution, &psi,
290: asls->dpsi, tao->stepdirection, &t, &ls_reason);
291: VecNorm(asls->dpsi, NORM_2, &ndpsi);
292: }
294: return(0);
295: }
297: /* ---------------------------------------------------------- */
301: PetscErrorCode TaoCreate_ASILS(TaoSolver tao)
302: {
303: TAO_SSLS *asls;
304: PetscErrorCode ierr;
305: const char *armijo_type = TAOLINESEARCH_ARMIJO;
308: PetscNewLog(tao,TAO_SSLS,&asls);
309: tao->data = (void*)asls;
310: tao->ops->solve = TaoSolve_ASILS;
311: tao->ops->setup = TaoSetUp_ASILS;
312: tao->ops->view = TaoView_SSLS;
313: tao->ops->setfromoptions = TaoSetFromOptions_SSLS;
314: tao->ops->destroy = TaoDestroy_ASILS;
315: tao->subset_type = TAO_SUBSET_SUBVEC;
316: asls->delta = 1e-10;
317: asls->rho = 2.1;
318: asls->fixed = PETSC_NULL;
319: asls->free = PETSC_NULL;
320: asls->J_sub = PETSC_NULL;
321: asls->Jpre_sub = PETSC_NULL;
322: asls->w = PETSC_NULL;
323: asls->r1 = PETSC_NULL;
324: asls->r2 = PETSC_NULL;
325: asls->r3 = PETSC_NULL;
326: asls->t1 = PETSC_NULL;
327: asls->t2 = PETSC_NULL;
328: asls->dxfree = PETSC_NULL;
329:
330: asls->identifier = 1e-5;
332: TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
333: TaoLineSearchSetType(tao->linesearch, armijo_type);
334: TaoLineSearchSetFromOptions(tao->linesearch);
336: KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
337: KSPSetFromOptions(tao->ksp);
338: tao->max_it = 2000;
339: tao->max_funcs = 4000;
340: tao->fatol = 0;
341: tao->frtol = 0;
342: tao->gttol = 0;
343: tao->grtol = 0;
344: tao->gatol = 1.0e-16;
345: tao->fmin = 1.0e-8;
347: return(0);
348: }