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: }