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 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(Tao tao) 59: {
60: TAO_SSLS *asls = (TAO_SSLS *)tao->data;
64: VecDuplicate(tao->solution,&tao->gradient);
65: VecDuplicate(tao->solution,&tao->stepdirection);
66: VecDuplicate(tao->solution,&asls->ff);
67: VecDuplicate(tao->solution,&asls->dpsi);
68: VecDuplicate(tao->solution,&asls->da);
69: VecDuplicate(tao->solution,&asls->db);
70: VecDuplicate(tao->solution,&asls->t1);
71: VecDuplicate(tao->solution,&asls->t2);
72: VecDuplicate(tao->solution, &asls->w);
73: asls->fixed = NULL;
74: asls->free = NULL;
75: asls->J_sub = NULL;
76: asls->Jpre_sub = NULL;
77: asls->r1 = NULL;
78: asls->r2 = NULL;
79: asls->r3 = NULL;
80: asls->dxfree = NULL;
81: return(0);
82: }
86: static PetscErrorCode Tao_ASLS_FunctionGradient(TaoLineSearch ls, Vec X, PetscReal *fcn, Vec G, void *ptr) 87: {
88: Tao tao = (Tao)ptr;
89: TAO_SSLS *asls = (TAO_SSLS *)tao->data;
93: TaoComputeConstraints(tao, X, tao->constraints);
94: VecFischer(X,tao->constraints,tao->XL,tao->XU,asls->ff);
95: VecNorm(asls->ff,NORM_2,&asls->merit);
96: *fcn = 0.5*asls->merit*asls->merit;
97: TaoComputeJacobian(tao,tao->solution,tao->jacobian,tao->jacobian_pre);
99: MatDFischer(tao->jacobian, tao->solution, tao->constraints,tao->XL, tao->XU, asls->t1, asls->t2,asls->da, asls->db);
100: VecPointwiseMult(asls->t1, asls->ff, asls->db);
101: MatMultTranspose(tao->jacobian,asls->t1,G);
102: VecPointwiseMult(asls->t1, asls->ff, asls->da);
103: VecAXPY(G,1.0,asls->t1);
104: return(0);
105: }
109: static PetscErrorCode TaoDestroy_ASFLS(Tao tao)110: {
111: TAO_SSLS *ssls = (TAO_SSLS *)tao->data;
115: VecDestroy(&ssls->ff);
116: VecDestroy(&ssls->dpsi);
117: VecDestroy(&ssls->da);
118: VecDestroy(&ssls->db);
119: VecDestroy(&ssls->w);
120: VecDestroy(&ssls->t1);
121: VecDestroy(&ssls->t2);
122: VecDestroy(&ssls->r1);
123: VecDestroy(&ssls->r2);
124: VecDestroy(&ssls->r3);
125: VecDestroy(&ssls->dxfree);
126: MatDestroy(&ssls->J_sub);
127: MatDestroy(&ssls->Jpre_sub);
128: ISDestroy(&ssls->fixed);
129: ISDestroy(&ssls->free);
130: PetscFree(tao->data);
131: tao->data = NULL;
132: return(0);
133: }
137: static PetscErrorCode TaoSolve_ASFLS(Tao tao)138: {
139: TAO_SSLS *asls = (TAO_SSLS *)tao->data;
140: PetscReal psi,ndpsi, normd, innerd, t=0;
141: PetscInt iter=0, nf;
142: PetscErrorCode ierr;
143: TaoConvergedReason reason;
144: TaoLineSearchConvergedReason ls_reason;
147: /* Assume that Setup has been called!
148: Set the structure for the Jacobian and create a linear solver. */
150: TaoComputeVariableBounds(tao);
151: TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,Tao_ASLS_FunctionGradient,tao);
152: TaoLineSearchSetObjectiveRoutine(tao->linesearch,Tao_SSLS_Function,tao);
153: TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);
155: VecMedian(tao->XL, tao->solution, tao->XU, tao->solution);
157: /* Calculate the function value and fischer function value at the
158: current iterate */
159: TaoLineSearchComputeObjectiveAndGradient(tao->linesearch,tao->solution,&psi,asls->dpsi);
160: VecNorm(asls->dpsi,NORM_2,&ndpsi);
162: while (1) {
163: /* Check the converged criteria */
164: PetscInfo3(tao,"iter %D, merit: %g, ||dpsi||: %g\n",iter, (double)asls->merit, (double)ndpsi);
165: TaoMonitor(tao, iter++, asls->merit, ndpsi, 0.0, t, &reason);
166: if (TAO_CONTINUE_ITERATING != reason) break;
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.
174: This rule gives superlinear asymptotic convergence
175: asls->atol = min(0.5, asls->merit*sqrt(asls->merit));
176: asls->rtol = 0.0;
178: This rule gives quadratic asymptotic convergence
179: asls->atol = min(0.5, asls->merit*asls->merit);
180: asls->rtol = 0.0;
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.
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: MatNorm(tao->jacobian,NORM_1,&asls->identifier);
191: asls->identifier = PetscMin(asls->merit, 1e-2) / (1 + asls->identifier);
193: VecSet(asls->t1,-asls->identifier);
194: VecSet(asls->t2, asls->identifier);
196: ISDestroy(&asls->fixed);
197: ISDestroy(&asls->free);
198: VecWhichBetweenOrEqual(asls->t1, asls->db, asls->t2, &asls->fixed);
199: ISComplementVec(asls->fixed,asls->t1, &asls->free);
201: ISGetSize(asls->fixed,&nf);
202: PetscInfo1(tao,"Number of fixed variables: %D\n", nf);
204: /* We now have our partition. Now calculate the direction in the
205: fixed variable space. */
206: TaoVecGetSubVec(asls->ff, asls->fixed, tao->subset_type, 0.0, &asls->r1);
207: TaoVecGetSubVec(asls->da, asls->fixed, tao->subset_type, 1.0, &asls->r2);
208: VecPointwiseDivide(asls->r1,asls->r1,asls->r2);
209: VecSet(tao->stepdirection,0.0);
210: VecISAXPY(tao->stepdirection, asls->fixed, 1.0,asls->r1);
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. */
217: TaoVecGetSubVec(asls->da, asls->free, tao->subset_type, 0.0, &asls->r1);
218: TaoVecGetSubVec(asls->ff, asls->free, tao->subset_type, 0.0, &asls->r2);
219: TaoVecGetSubVec(asls->db, asls->free, tao->subset_type, 1.0, &asls->r3);
220: VecPointwiseDivide(asls->r1,asls->r1, asls->r3);
221: VecPointwiseDivide(asls->r2,asls->r2, asls->r3);
223: /* r1 is the diagonal perturbation
224: r2 is the right hand side
225: r3 is no longer needed
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. */
231: MatMult(tao->jacobian, tao->stepdirection, asls->t1);
232: TaoVecGetSubVec(asls->t1,asls->free,tao->subset_type,0.0,&asls->r3);
233: VecAXPY(asls->r2, -1.0, asls->r3);
235: /* Calculate the reduced problem matrix and the direction */
236: TaoMatGetSubMat(tao->jacobian, asls->free, asls->w, tao->subset_type,&asls->J_sub);
237: if (tao->jacobian != tao->jacobian_pre) {
238: TaoMatGetSubMat(tao->jacobian_pre, asls->free, asls->w, tao->subset_type, &asls->Jpre_sub);
239: } else {
240: MatDestroy(&asls->Jpre_sub);
241: asls->Jpre_sub = asls->J_sub;
242: PetscObjectReference((PetscObject)(asls->Jpre_sub));
243: }
244: MatDiagonalSet(asls->J_sub, asls->r1,ADD_VALUES);
245: TaoVecGetSubVec(tao->stepdirection, asls->free, tao->subset_type, 0.0, &asls->dxfree);
246: VecSet(asls->dxfree, 0.0);
248: /* Calculate the reduced direction. (Really negative of Newton
249: direction. Therefore, rest of the code uses -d.) */
250: KSPReset(tao->ksp);
251: KSPSetOperators(tao->ksp, asls->J_sub, asls->Jpre_sub);
252: KSPSolve(tao->ksp, asls->r2, asls->dxfree);
254: /* Add the direction in the free variables back into the real direction. */
255: VecISAXPY(tao->stepdirection, asls->free, 1.0,asls->dxfree);
258: /* Check the projected real direction for descent and if not, use the negative
259: gradient direction. */
260: VecCopy(tao->stepdirection, asls->w);
261: VecScale(asls->w, -1.0);
262: VecBoundGradientProjection(asls->w, tao->solution, tao->XL, tao->XU, asls->w);
263: VecNorm(asls->w, NORM_2, &normd);
264: VecDot(asls->w, asls->dpsi, &innerd);
266: if (innerd >= -asls->delta*pow(normd, asls->rho)) {
267: PetscInfo1(tao,"Gradient direction: %5.4e.\n", (double)innerd);
268: PetscInfo1(tao, "Iteration %D: newton direction not descent\n", iter);
269: VecCopy(asls->dpsi, tao->stepdirection);
270: VecDot(asls->dpsi, tao->stepdirection, &innerd);
271: }
273: VecScale(tao->stepdirection, -1.0);
274: innerd = -innerd;
276: /* We now have a correct descent direction. Apply a linesearch to
277: find the new iterate. */
278: TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0);
279: TaoLineSearchApply(tao->linesearch, tao->solution, &psi,
280: asls->dpsi, tao->stepdirection, &t, &ls_reason);
281: VecNorm(asls->dpsi, NORM_2, &ndpsi);
282: }
283: return(0);
284: }
286: /* ---------------------------------------------------------- */
287: /*MC
288: TAOASFLS - Active-set feasible linesearch algorithm for solving
289: complementarity constraints
291: Options Database Keys:
292: + -tao_ssls_delta - descent test fraction
293: - -tao_ssls_rho - descent test power
295: Level: beginner
296: M*/
297: EXTERN_C_BEGIN
300: PetscErrorCode TaoCreate_ASFLS(Tao tao)301: {
302: TAO_SSLS *asls;
304: const char *armijo_type = TAOLINESEARCHARMIJO;
307: PetscNewLog(tao,&asls);
308: tao->data = (void*)asls;
309: tao->ops->solve = TaoSolve_ASFLS;
310: tao->ops->setup = TaoSetUp_ASFLS;
311: tao->ops->view = TaoView_SSLS;
312: tao->ops->setfromoptions = TaoSetFromOptions_SSLS;
313: tao->ops->destroy = TaoDestroy_ASFLS;
314: tao->subset_type = TAO_SUBSET_SUBVEC;
315: asls->delta = 1e-10;
316: asls->rho = 2.1;
317: asls->fixed = NULL;
318: asls->free = NULL;
319: asls->J_sub = NULL;
320: asls->Jpre_sub = NULL;
321: asls->w = NULL;
322: asls->r1 = NULL;
323: asls->r2 = NULL;
324: asls->r3 = NULL;
325: asls->t1 = NULL;
326: asls->t2 = NULL;
327: asls->dxfree = NULL;
328: asls->identifier = 1e-5;
330: TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
331: TaoLineSearchSetType(tao->linesearch, armijo_type);
332: TaoLineSearchSetFromOptions(tao->linesearch);
334: KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
335: KSPSetFromOptions(tao->ksp);
336: tao->max_it = 2000;
337: tao->max_funcs = 4000;
338: tao->fatol = 0;
339: tao->frtol = 0;
340: tao->gttol = 0;
341: tao->grtol = 0;
342: #if defined(PETSC_USE_REAL_SINGLE)
343: tao->gatol = 1.0e-6;
344: tao->fmin = 1.0e-4;
345: #else
346: tao->gatol = 1.0e-16;
347: tao->fmin = 1.0e-8;
348: #endif
351: return(0);
352: }
353: EXTERN_C_END