Actual source code: asfls.c

petsc-3.6.1 2015-08-06
Report Typos and Errors
  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                     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",tao->niter,(double)asls->merit,(double)ndpsi);
165:     TaoMonitor(tao,tao->niter,asls->merit,ndpsi,0.0,t,&reason);
166:     if (TAO_CONTINUE_ITERATING != reason) break;
167:     tao->niter++;

169:     /* We are going to solve a linear system of equations.  We need to
170:        set the tolerances for the solve so that we maintain an asymptotic
171:        rate of convergence that is superlinear.
172:        Note: these tolerances are for the reduced system.  We really need
173:        to make sure that the full system satisfies the full-space conditions.

175:        This rule gives superlinear asymptotic convergence
176:        asls->atol = min(0.5, asls->merit*sqrt(asls->merit));
177:        asls->rtol = 0.0;

179:        This rule gives quadratic asymptotic convergence
180:        asls->atol = min(0.5, asls->merit*asls->merit);
181:        asls->rtol = 0.0;

183:        Calculate a free and fixed set of variables.  The fixed set of
184:        variables are those for the d_b is approximately equal to zero.
185:        The definition of approximately changes as we approach the solution
186:        to the problem.

188:        No one rule is guaranteed to work in all cases.  The following
189:        definition is based on the norm of the Jacobian matrix.  If the
190:        norm is large, the tolerance becomes smaller. */
191:     MatNorm(tao->jacobian,NORM_1,&asls->identifier);
192:     asls->identifier = PetscMin(asls->merit, 1e-2) / (1 + asls->identifier);

194:     VecSet(asls->t1,-asls->identifier);
195:     VecSet(asls->t2, asls->identifier);

197:     ISDestroy(&asls->fixed);
198:     ISDestroy(&asls->free);
199:     VecWhichBetweenOrEqual(asls->t1, asls->db, asls->t2, &asls->fixed);
200:     ISComplementVec(asls->fixed,asls->t1, &asls->free);

202:     ISGetSize(asls->fixed,&nf);
203:     PetscInfo1(tao,"Number of fixed variables: %D\n", nf);

205:     /* We now have our partition.  Now calculate the direction in the
206:        fixed variable space. */
207:     TaoVecGetSubVec(asls->ff, asls->fixed, tao->subset_type, 0.0, &asls->r1);
208:     TaoVecGetSubVec(asls->da, asls->fixed, tao->subset_type, 1.0, &asls->r2);
209:     VecPointwiseDivide(asls->r1,asls->r1,asls->r2);
210:     VecSet(tao->stepdirection,0.0);
211:     VecISAXPY(tao->stepdirection, asls->fixed, 1.0,asls->r1);

213:     /* Our direction in the Fixed Variable Set is fixed.  Calculate the
214:        information needed for the step in the Free Variable Set.  To
215:        do this, we need to know the diagonal perturbation and the
216:        right hand side. */

218:     TaoVecGetSubVec(asls->da, asls->free, tao->subset_type, 0.0, &asls->r1);
219:     TaoVecGetSubVec(asls->ff, asls->free, tao->subset_type, 0.0, &asls->r2);
220:     TaoVecGetSubVec(asls->db, asls->free, tao->subset_type, 1.0, &asls->r3);
221:     VecPointwiseDivide(asls->r1,asls->r1, asls->r3);
222:     VecPointwiseDivide(asls->r2,asls->r2, asls->r3);

224:     /* r1 is the diagonal perturbation
225:        r2 is the right hand side
226:        r3 is no longer needed

228:        Now need to modify r2 for our direction choice in the fixed
229:        variable set:  calculate t1 = J*d, take the reduced vector
230:        of t1 and modify r2. */

232:     MatMult(tao->jacobian, tao->stepdirection, asls->t1);
233:     TaoVecGetSubVec(asls->t1,asls->free,tao->subset_type,0.0,&asls->r3);
234:     VecAXPY(asls->r2, -1.0, asls->r3);

236:     /* Calculate the reduced problem matrix and the direction */
237:     TaoMatGetSubMat(tao->jacobian, asls->free, asls->w, tao->subset_type,&asls->J_sub);
238:     if (tao->jacobian != tao->jacobian_pre) {
239:       TaoMatGetSubMat(tao->jacobian_pre, asls->free, asls->w, tao->subset_type, &asls->Jpre_sub);
240:     } else {
241:       MatDestroy(&asls->Jpre_sub);
242:       asls->Jpre_sub = asls->J_sub;
243:       PetscObjectReference((PetscObject)(asls->Jpre_sub));
244:     }
245:     MatDiagonalSet(asls->J_sub, asls->r1,ADD_VALUES);
246:     TaoVecGetSubVec(tao->stepdirection, asls->free, tao->subset_type, 0.0, &asls->dxfree);
247:     VecSet(asls->dxfree, 0.0);

249:     /* Calculate the reduced direction.  (Really negative of Newton
250:        direction.  Therefore, rest of the code uses -d.) */
251:     KSPReset(tao->ksp);
252:     KSPSetOperators(tao->ksp, asls->J_sub, asls->Jpre_sub);
253:     KSPSolve(tao->ksp, asls->r2, asls->dxfree);
254:     KSPGetIterationNumber(tao->ksp,&tao->ksp_its);
255:     tao->ksp_tot_its+=tao->ksp_its;

257:     /* Add the direction in the free variables back into the real direction. */
258:     VecISAXPY(tao->stepdirection, asls->free, 1.0,asls->dxfree);


261:     /* Check the projected real direction for descent and if not, use the negative
262:        gradient direction. */
263:     VecCopy(tao->stepdirection, asls->w);
264:     VecScale(asls->w, -1.0);
265:     VecBoundGradientProjection(asls->w, tao->solution, tao->XL, tao->XU, asls->w);
266:     VecNorm(asls->w, NORM_2, &normd);
267:     VecDot(asls->w, asls->dpsi, &innerd);

269:     if (innerd >= -asls->delta*PetscPowReal(normd, asls->rho)) {
270:       PetscInfo1(tao,"Gradient direction: %5.4e.\n", (double)innerd);
271:       PetscInfo1(tao, "Iteration %D: newton direction not descent\n", tao->niter);
272:       VecCopy(asls->dpsi, tao->stepdirection);
273:       VecDot(asls->dpsi, tao->stepdirection, &innerd);
274:     }

276:     VecScale(tao->stepdirection, -1.0);
277:     innerd = -innerd;

279:     /* We now have a correct descent direction.  Apply a linesearch to
280:        find the new iterate. */
281:     TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0);
282:     TaoLineSearchApply(tao->linesearch, tao->solution, &psi,asls->dpsi, tao->stepdirection, &t, &ls_reason);
283:     VecNorm(asls->dpsi, NORM_2, &ndpsi);
284:   }
285:   return(0);
286: }

288: /* ---------------------------------------------------------- */
289: /*MC
290:    TAOASFLS - Active-set feasible linesearch algorithm for solving
291:        complementarity constraints

293:    Options Database Keys:
294: + -tao_ssls_delta - descent test fraction
295: - -tao_ssls_rho - descent test power

297:    Level: beginner
298: M*/
301: PETSC_EXTERN PetscErrorCode TaoCreate_ASFLS(Tao tao)
302: {
303:   TAO_SSLS       *asls;
305:   const char     *armijo_type = TAOLINESEARCHARMIJO;

308:   PetscNewLog(tao,&asls);
309:   tao->data = (void*)asls;
310:   tao->ops->solve = TaoSolve_ASFLS;
311:   tao->ops->setup = TaoSetUp_ASFLS;
312:   tao->ops->view = TaoView_SSLS;
313:   tao->ops->setfromoptions = TaoSetFromOptions_SSLS;
314:   tao->ops->destroy = TaoDestroy_ASFLS;
315:   tao->subset_type = TAO_SUBSET_SUBVEC;
316:   asls->delta = 1e-10;
317:   asls->rho = 2.1;
318:   asls->fixed = NULL;
319:   asls->free = NULL;
320:   asls->J_sub = NULL;
321:   asls->Jpre_sub = NULL;
322:   asls->w = NULL;
323:   asls->r1 = NULL;
324:   asls->r2 = NULL;
325:   asls->r3 = NULL;
326:   asls->t1 = NULL;
327:   asls->t2 = NULL;
328:   asls->dxfree = NULL;
329:   asls->identifier = 1e-5;

331:   TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
332:   TaoLineSearchSetType(tao->linesearch, armijo_type);
333:   TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
334:   TaoLineSearchSetFromOptions(tao->linesearch);

336:   KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
337:   KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix);
338:   KSPSetFromOptions(tao->ksp);

340:   /* Override default settings (unless already changed) */
341:   if (!tao->max_it_changed) tao->max_it = 2000;
342:   if (!tao->max_funcs_changed) tao->max_funcs = 4000;
343:   if (!tao->fatol_changed) tao->fatol = 0;
344:   if (!tao->frtol_changed) tao->frtol = 0;
345:   if (!tao->gttol_changed) tao->gttol = 0;
346:   if (!tao->grtol_changed) tao->grtol = 0;
347: #if defined(PETSC_USE_REAL_SINGLE)
348:   if (!tao->gatol_changed) tao->gatol = 1.0e-6;
349:   if (!tao->fmin_changed)  tao->fmin = 1.0e-4;
350: #else
351:   if (!tao->gatol_changed) tao->gatol = 1.0e-16;
352:   if (!tao->fmin_changed)  tao->fmin = 1.0e-8;
353: #endif
354:   return(0);
355: }