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