Actual source code: isils.c

      1: #include "src/complementarity/impls/ssls/ssls.h"

      3: // Context for ISXLS
      4: //   -- interior-point  - smoothed fisched function
      5: //   -- active-set      - reduced matrices formed
      6: //                          - inherit properties of original system
      7: //   -- semismooth (S)  - function not differentiable
      8: //                      - merit function continuously differentiable
      9: //                      - Fischer-Burmeister reformulation of complementarity
     10: //                        - Billups composition for two finite bounds
     11: //   -- infeasible (I)  - iterates not guaranteed to remain within bounds
     12: //   -- feasible (F)    - iterates guaranteed to remain within bounds
     13: //   -- linesearch (LS) - Armijo rule on direction
     14: //
     15: // Many other reformulations are possible and combinations of
     16: // feasible/infeasible and linesearch/trust region are possible.
     17: //
     18: // Basic theory
     19: //   Fischer-Burmeister reformulation is semismooth with a continuously
     20: //   differentiable merit function and strongly semismooth if the F has
     21: //   lipschitz continuous derivatives.
     22: //
     23: //   Every accumulation point generated by the algorithm is a stationary
     24: //   point for the merit function.  Stationary points of the merit function
     25: //   are solutions of the complementarity problem if
     26: //     a.  the stationary point has a BD-regular subdifferential, or
     27: //     b.  the Schur complement F'/F'_ff is a P_0-matrix where ff is the
     28: //         index set corresponding to the free variables.
     29: //
     30: //   If one of the accumulation points has a BD-regular subdifferential then
     31: //     a.  the entire sequence converges to this accumulation point at
     32: //         a local q-superlinear rate
     33: //     b.  if in addition the reformulation is strongly semismooth near
     34: //         this accumulation point, then the algorithm converges at a
     35: //         local q-quadratic rate.
     36: //
     37: // The theory for the feasible version follows from the feasible descent
     38: // algorithm framework.
     39: //
     40: // References:
     41: //   Billups, "Algorithms for Complementarity Problems and Generalized
     42: //     Equations," Ph.D thesis, University of Wisconsin - Madison, 1995.
     43: //   De Luca, Facchinei, Kanzow, "A Semismooth Equation Approach to the
     44: //     Solution of Nonlinear Complementarity Problems," Mathematical
     45: //     Programming, 75, pages 407-439, 1996.
     46: //   Ferris, Kanzow, Munson, "Feasible Descent Algorithms for Mixed
     47: //     Complementarity Problems," Mathematical Programming, 86,
     48: //     pages 475-497, 1999.
     49: //   Fischer, "A Special Newton-type Optimization Method," Optimization,
     50: //     24, pages 269-284, 1992
     51: //   Munson, Facchinei, Ferris, Fischer, Kanzow, "The Semismooth Algorithm
     52: //     for Large Scale Complementarity Problems," Technical Report 99-06,
     53: //     University of Wisconsin - Madison, 1999.

     55: /*------------------------------------------------------------*/
     58: int Tao_ISLS_SFunction(TAO_SOLVER tao, TaoVec *X, double *fcn, void *solver)
     59: {
     60:   TAO_SSLS *isls = (TAO_SSLS *)solver;
     61:   TaoVec *l, *u;
     62:   int info;

     64:   info = TaoGetVariableBounds(tao, &l, &u); CHKERRQ(info);

     66:   info = TaoComputeConstraints(tao, X, isls->f); CHKERRQ(info);
     67:   info = isls->ff->SFischer(X, isls->f, l, u, isls->mu); CHKERRQ(info);
     68:   info = isls->ff->Norm2squared(&isls->merit_eqn); CHKERRQ(info);

     70:   // Add in the constraint for mu
     71:   isls->mucon = exp(isls->mu) - 1.0;
     72:   isls->merit_mu = isls->mucon*isls->mucon;

     74:   isls->merit = isls->merit_eqn + isls->merit_mu;

     76:   // Calculate the objective and merit function
     77:   *fcn = 0.5*isls->merit;
     78:   isls->merit = sqrt(isls->merit);
     79:   return 0;
     80: }

     82: /*------------------------------------------------------------*/
     85: int Tao_ISLS_SFunctionGradient(TAO_SOLVER tao, TaoVec *X, double *fcn,
     86:                                TaoVec *G, void *solver)
     87: {
     88:   TAO_SSLS *isls = (TAO_SSLS *)solver;
     89:   TaoVec *l, *u;
     90:   TaoMat *J;
     91:   int info;

     93:   info = TaoGetVariableBounds(tao, &l, &u); CHKERRQ(info);
     94:   info = TaoGetJacobian(tao, &J); CHKERRQ(info);

     96:   info = TaoComputeConstraints(tao, X, isls->f); CHKERRQ(info);
     97:   info = isls->ff->SFischer(X, isls->f, l, u, isls->mu); CHKERRQ(info);
     98:   info = isls->ff->Norm2squared(&isls->merit_eqn); CHKERRQ(info);

    100:   // Add in the constraint for mu
    101:   isls->mucon = exp(isls->mu) - 1.0;
    102:   isls->merit_mu = isls->mucon*isls->mucon;

    104:   isls->merit = isls->merit_eqn + isls->merit_mu;

    106:   // Calculate the objective and merit function
    107:   *fcn = 0.5*isls->merit;
    108:   isls->merit = sqrt(isls->merit);

    110:   info = TaoComputeJacobian(tao, X, J); CHKERRQ(info);
    111:   info = J->D_SFischer(X, isls->f, l, u, isls->mu,
    112:                        isls->t1, isls->t2,
    113:                        isls->da, isls->db, isls->dm); CHKERRQ(info);
    114:   isls->d_mucon = exp(isls->mu);

    116:   info = isls->t1->PointwiseMultiply(isls->ff, isls->db);
    117:   info = J->MultiplyTranspose(isls->t1, G); CHKERRQ(info);
    118:   info = isls->t1->PointwiseMultiply(isls->ff, isls->da); CHKERRQ(info);
    119:   info = G->Axpy(1.0, isls->t1); CHKERRQ(info);
    120:   info = isls->dm->Dot(isls->ff, &isls->g_mucon); CHKERRQ(info);
    121:   isls->g_mucon += exp(isls->mu)*isls->mucon;
    122:   return 0;
    123: }

    125: /*------------------------------------------------------------*/
    128: static int TaoSolve_ISILS(TAO_SOLVER tao, void *solver)
    129: {
    130:   TAO_SSLS *isls = (TAO_SSLS *)solver;

    132:   TaoVec *x;                        // Solution vector
    133:   TaoVec *l, *u;                // Problem bounds
    134:   TaoMat *J;                        // Jacobian matrix

    136:   TaoIndexSet *FreeVariableSet;        // Free variables
    137:   TaoIndexSet *FixedVariableSet;// Fixed variables

    139:   TaoMat *Jsub;                        // Reduced Jacobian
    140:   TaoVec *r1, *r2, *r3;                // Temporary vectors

    142: //  TaoLinearSolver *lsolver;

    144:   double psi, ndpsi, normd, innerd, t=0;
    145:   double two_norm;
    146:   double sigma = 1e-4, beta = 0.5, minstep = TAO_EPSILON;
    147:   double muval, psival;
    148:   int iter=0, step_type, info, nn, nf;
    149:   TaoTerminateReason reason;
    150:   TaoTruth flag;

    152:   TaoFunctionBegin;

    154:   // Assume that Setup has been called!

    156:   // Start by obtaining the solution vector, Jacobian,
    157:   // and problem size
    158:   info = TaoGetSolution(tao, &x); CHKERRQ(info);
    159:   info = TaoGetJacobian(tao, &J); CHKERRQ(info);
    160:   info = x->GetDimension(&nn); CHKERRQ(info);

    162:   // Evaluate the bounds
    163:   info = TaoGetVariableBounds(tao, &l, &u); CHKERRQ(info);
    164:   info = TaoEvaluateVariableBounds(tao, l, u); CHKERRQ(info);

    166:   // Create vector for derivative wrt smoothing parameter
    167:   info = x->Clone(&(isls->dm)); CHKERRQ(info);

    169:   // Create temporary vectors and index sets (for active set)
    170:   info = x->CreateIndexSet(&FreeVariableSet); CHKERRQ(info);
    171:   info = x->CreateIndexSet(&FixedVariableSet); CHKERRQ(info);
    172:   info = x->Clone(&r1); CHKERRQ(info);
    173:   info = x->Clone(&r2); CHKERRQ(info);
    174:   info = x->Clone(&r3); CHKERRQ(info);

    176:   // Create submatrix
    177:   info = J->CreateReducedMatrix(FreeVariableSet, FreeVariableSet, &Jsub); CHKERRQ(info);

    179:   // Initialize the merit function
    180:   //   Used in the linesearch
    181:   //   Termination based on different function
    182:   info = TaoSetMeritFunction(tao, Tao_ISLS_SFunction, Tao_ISLS_SFunctionGradient,
    183:                              TAO_NULL, TAO_NULL, isls); CHKERRQ(info);

    185:   // Initialize the value for the smoothing parameter
    186:   isls->mu = isls->mu_init;

    188: #if 0
    189:   // Choice of mu:
    190:   //   Simple search

    192:   info = TaoComputeMeritFunctionGradient(tao, x, &psi, isls->dpsi); CHKERRQ(info);

    194:   psival = fabs(isls->merit_eqn - nn*isls->merit_mu);
    195:   printf("mu: %5.4e merit_eqn: %5.4e merit_mu: %5.4e\n",
    196:          isls->mu, isls->merit_eqn, nn*isls->merit_mu);

    198:   double fact;
    199:   if (isls->merit_eqn < nn*isls->merit_mu) {
    200:     fact = 0.5;
    201:   } else
    202:     {
    203:       fact = 2.0;
    204:     }

    206:   isls->mu *= fact;
    207:   info = TaoComputeMeritFunctionGradient(tao, x, &psi, isls->dpsi); CHKERRQ(info);
    208:   printf("mu: %5.4e merit_eqn: %5.4e merit_mu: %5.4e\n",
    209:          isls->mu, isls->merit_eqn, nn*isls->merit_mu);
    210:
    211:   while (fabs(isls->merit_eqn - nn*isls->merit_mu) < psival) {
    212:     psival = fabs(isls->merit_eqn - nn*isls->merit_mu);

    214:     isls->mu *= fact;
    215:     info = TaoComputeMeritFunctionGradient(tao, x, &psi, isls->dpsi); CHKERRQ(info);
    216:     printf("mu: %5.4e merit_eqn: %5.4e merit_mu: %5.4e\n",
    217:            isls->mu, isls->merit_eqn, nn*isls->merit_mu);
    218:   }

    220:   isls->mu /= fact;
    221:   printf("mu_choice : %5.4e\n", isls->mu);
    222: #endif

    224:   // DONE:
    225:   //   1.  Eliminated smoothing parameter vector.  Do not know how to deal
    226:   //       with mixture of differentiable and nondifferentiable functions.
    227:   //       Maybe let the smoothing parameter be a vector later.
    228:   //   2.  Rewrote smoothed functions to provide the following information:
    229:   //       - derivative of the function with respect to smoothing parameter
    230:   //       - gives nonsmooth function when smoothing parameter == 0
    231:   //   3.  Rewrote the function and function/gradient routines to add
    232:   //       the smoothing constraint (e^{mu}-1 = 0) and the appropriate
    233:   //       derivatives.
    234:   //   4.  Modified the right-hand side of the direction calculation to
    235:   //       account for the smoothing constraint.  I.e. the direction in
    236:   //       the smoothing parameter is substituted out of the problem.
    237:   //   5.  Modified descent direction test to deal with the augmented
    238:   //       problem with the smoothing parameter.
    239:   //   6.  Hacked in a monotone linesearch to deal with the augmented system.
    240:   //   7.  When smoothing parameter == 0, make sure we get the same
    241:   //       behavior as asils.
    242:   // TODO:
    243:   //   8.  Add method to optimally choose initial value for mu.

    245:
    246:   // Calculate the fischer function at the initial iterate
    247:   Tao_ASLS_FunctionGradient(tao, x, &psi, isls->dpsi, isls);
    248:   info = isls->dpsi->Norm2(&ndpsi); CHKERRQ(info);

    250:   while (1) {

    252:     // Assert:
    253:     //   Fischer function for terminatation
    254:     //   Evaluated at current iterate
    255:     //   Norm gradient of Fischer function calculated

    257:     printf("nd_merit : %5.4e\n", isls->merit);
    258:     printf("nd_ndpsi : %5.4e\n", ndpsi);

    260:     // Check the termination criteria

    262:     info = TaoMonitor(tao, iter++, isls->merit, ndpsi, 0.0, t, &reason); CHKERRQ(info);
    263:     if (TAO_CONTINUE_ITERATING != reason) break;

    265:     // Compute the smoothed Fischer function.
    266:     // This forms the basis for the interior-point method.

    268:     info = TaoComputeMeritFunctionGradient(tao, x, &psi, isls->dpsi); CHKERRQ(info);

    270:     printf("merit    : %5.4e\n", isls->merit);
    271:     printf("merit_eqn: %5.4e\n", isls->merit_eqn);
    272:     printf("merit_mu : %5.4e\n", isls->merit_mu);
    273:     printf("mu       : %5.4e\n", isls->mu);
    274:     printf("mucon    : %5.4e\n", isls->mucon);
    275:     printf("d_mucon  : %5.4e\n", isls->d_mucon);
    276:     printf("g_mucon  : %5.4e\n", isls->g_mucon);

    278: #if 0
    279:     // Compute tolerances for linear system solver

    281:     // We want to set tolerances so that we maintain a superlinear
    282:     // asymptotic rate of convergence.  Note: the tolerances set are
    283:     // for the reduced system.  We really need to make sure that the
    284:     // full system satisfies the termination conditions.

    286:     // Superlinear asymptotic convergence:
    287:     // isls->atol = TaoMin(0.5, isls->merit*sqrt(isls->merit));
    288:     // isls->rtol = 0.0;

    290:     // Quadratic asymptotic convergence:
    291:     isls->atol = TaoMin(0.5, isls->merit*isls->merit);
    292:     isls->rtol = 0.0;

    294:     // Set the convergence tolerances.
    295:     info = TaoGetLinearSolver(tao, &lsolver); CHKERRQ(info);
    296:     // info = lsolver->SetTolerances(isls->rtol, isls->atol, 1e30, 10*nn); CHKERRQ(info);
    297:     // info = lsolver->SetTolerances(0.0, 0.0, 1e30, 10*nn); CHKERRQ(info);
    298:     // info = lsolver->SetTolerances(0.0, 0.0, 1e30, 10*nn); CHKERRQ(info);
    299: #endif

    301:     // Determine whether we are taking a predictor or corrector step
    302:     // Calculate the direction in the smoothing parameter
    303:     // Substitute out of the problem

    305:     if ((isls->merit_eqn <= 1e-10) ||
    306:         (isls->merit_eqn <= nn*isls->merit_mu)) {
    307:       // Only take a predictor step when ||eqn|| <= const*||mu||
    308:       printf("Predictor step.\n");
    309:       step_type = 1;
    310:       isls->dmu = isls->mucon / isls->d_mucon;
    311:       printf("dmu    : %5.4e\n", isls->dmu);
    312:     }
    313:     else {
    314:       // Otherwise take a correct step
    315:       printf("Corrector step.\n");
    316:       step_type = 0;
    317:       isls->dmu = 0;
    318:     }
    319:     info = isls->w->Waxpby(1.0, isls->ff, -isls->dmu, isls->dm);

    321:     // Assert: w vector now contains right-hand side
    322:     // We solve the system:
    323:     //    [DPhi DMu ][dx ] = Phi
    324:     //    [   0 DRho][dmu] = Rho
    325:     // where Rho is the forcing function.

    327:     // Now calculate a free and fixed set of variables.  The fixed set of
    328:     // variables are those for which d_b is approximately equal to zero.
    329:     // The definition of approximately changes as we approact the
    330:     // solution.

    332:     // No one rule is guaranteed to work in all cases.  The following
    333:     // definition is based on the norm of the Jacobian matrix.  If the
    334:     // norm is large, the tolerance becomes smaller.

    336:     info = J->Norm1(&(isls->identifier)); CHKERRQ(info);
    337:     isls->identifier = TaoMin(isls->merit, 1e-2) / (1 + isls->identifier);

    339:     info = isls->t1->SetToConstant(-isls->identifier); CHKERRQ(info);
    340:     info = isls->t2->SetToConstant( isls->identifier); CHKERRQ(info);

    342:     info = isls->t1->SetToConstant(0.0); CHKERRQ(info);
    343:     info = isls->t2->SetToConstant(0.0); CHKERRQ(info);

    345:     info = FixedVariableSet->WhichBetweenOrEqual(isls->t1, isls->db, isls->t2); CHKERRQ(info);
    346:     info = FreeVariableSet->ComplementOf(FixedVariableSet); CHKERRQ(info);

    348:     info = FixedVariableSet->GetSize(&nf);
    349:     if (nf) printf("Fixed size: %d\n", nf);

    351:     // Assert:
    352:     //   Partition created.
    353:     // Now calculate the direction in the fixed variable space.
    354:     info = r1->SetReducedVec(isls->w, FixedVariableSet); CHKERRQ(info);
    355:     info = r2->SetReducedVec(isls->da, FixedVariableSet); CHKERRQ(info);
    356:     info = r1->PointwiseDivide(r1, r2); CHKERRQ(info);

    358:     // Throw the fixed variable direction into a global direction
    359:     info = isls->d->SetToZero(); CHKERRQ(info);
    360:     info = isls->d->ReducedXPY(r1, FixedVariableSet); CHKERRQ(info);

    362:     // Our direction in the FixedVariableSet is fixed.  Calculate the
    363:     // information needed for the step in the FreeVariableSet.  To
    364:     // do this, we need to know the diagonal perturbation and the
    365:     // right hand side.
    366:     info = r1->SetReducedVec(isls->da, FreeVariableSet); CHKERRQ(info);
    367:     info = r2->SetReducedVec(isls->w,  FreeVariableSet); CHKERRQ(info);
    368:     info = r3->SetReducedVec(isls->db, FreeVariableSet); CHKERRQ(info);

    370:     info = r1->PointwiseDivide(r1, r3); CHKERRQ(info);
    371:     info = r2->PointwiseDivide(r2, r3); CHKERRQ(info);

    373:     // Assert:
    374:     //   r1 is the diagonal perturbation
    375:     //   r2 is the right hand side
    376:     //   r3 is no longer needed
    377:     // Now need to modify r2 for our direction choice in the fixed
    378:     // variable set:  calculate t1 = J*d, take the reduced vector
    379:     // of t1 and modify r2.
    380:     info = J->Multiply(isls->d, isls->t1); CHKERRQ(info);
    381:     info = r3->SetReducedVec(isls->t1, FreeVariableSet); CHKERRQ(info);
    382:     info = r2->Axpy(-1.0, r3);

    384:     // Calculate the reduced problem matrix
    385:     info = Jsub->SetReducedMatrix(J, FreeVariableSet, FreeVariableSet); CHKERRQ(info);
    386:     info = Jsub->AddDiagonal(r1); CHKERRQ(info);

    388:     // Assert:
    389:     //   r1 is no longer needed
    390:     // Calculate the reduced direction.  (Really negative of Newton
    391:     // direction.  Therefore, rest of the code uses -d.)
    392:     info = TaoPreLinearSolve(tao, Jsub); CHKERRQ(info);
    393:     info = TaoLinearSolve(tao, Jsub, r2, r1, &flag); CHKERRQ(info);

    395:     // Add the direction in the free variables back into the real direction.
    396:     info = isls->d->ReducedXPY(r1, FreeVariableSet); CHKERRQ(info);

    398:     // Calculate the norm of the relative residual in the real space
    399:     info = isls->t1->PointwiseMultiply(isls->da, isls->d); CHKERRQ(info);
    400:     info = J->Multiply(isls->d, isls->t2); CHKERRQ(info);
    401:     info = isls->t2->PointwiseMultiply(isls->db, isls->t2); CHKERRQ(info);
    402:     info = isls->t1->Axpy( 1.0, isls->t2); CHKERRQ(info);
    403:     info = isls->t1->Axpy(-1.0, isls->w); CHKERRQ(info);
    404:     info = isls->t1->Norm2(&two_norm); CHKERRQ(info);

    406:     printf("eqn_resid: %5.4e\n", two_norm);

    408:     // Check direction for descent.
    409:     info = isls->d->Norm2squared(&normd); CHKERRQ(info);
    410:     normd += isls->dmu*isls->dmu;
    411:     normd = sqrt(normd);

    413:     info = isls->d->Dot(isls->dpsi, &innerd); CHKERRQ(info);
    414:     innerd += isls->g_mucon*isls->dmu;

    416:     printf("normd  : %5.4e\n", normd);
    417:     printf("innerd : %5.4e\n", innerd);

    419:     if (innerd <= isls->delta*pow(normd, isls->rho)) {
    420:       // Descent direction test failed, use gradient direction
    421:       info = PetscLogInfo((tao, "TaoSolve_ISILS: %d: newton direction not descent\n", iter)); CHKERRQ(info);
    422:       printf("Gradient direction\n");

    424:       info = isls->d->CopyFrom(isls->dpsi); CHKERRQ(info);
    425:       if (1 == step_type) {
    426:         // Predictor step gradient
    427:         isls->dmu = isls->g_mucon;
    428:       }
    429:       else {
    430:         // Corrector step gradient
    431:         isls->dmu = 0;
    432:       }

    434:       info = isls->d->Dot(isls->dpsi, &innerd); CHKERRQ(info);
    435:       innerd += isls->g_mucon*isls->dmu;

    437:       printf("g_innerd : %5.4e\n", innerd);
    438:     }

    440:     // We have a descent direction.  We need to take the negative to
    441:     // correct the sign of the direction.
    442:     info = isls->d->Negate(); CHKERRQ(info);
    443:     isls->dmu = -isls->dmu;
    444:     innerd = -innerd;

    446:     // Apply an armijo linesearch to find the new iterate.
    447:     innerd *= sigma;
    448:     muval = isls->mu;
    449:     t = 1;

    451:     while (t >= minstep) {
    452:       // Calculate iterate
    453:       info = isls->w->Waxpby(1.0, x, t, isls->d); CHKERRQ(info);
    454:       isls->mu = muval + t*isls->dmu;
    455:
    456:       // Calculate function at new iterate
    457:       info = TaoComputeMeritFunction(tao, isls->w, &psival); CHKERRQ(info);

    459:       printf("psival: %5.4e psi+t*innerd: %5.4e\n", psival, psi + t*innerd);

    461:       // Check descent condition
    462:       if (psival <= psi + t*innerd) {
    463:         break;
    464:       }
    465:       t *= beta;
    466:     }

    468:     x->CopyFrom(isls->w);

    470:     printf("t      : %5.4e\n", t);
    471:     printf("psi    : %5.4e\n", psi);
    472:     printf("psival : %5.4e\n", psival);
    473:     printf("diff   : %5.4e\n", psival - psi);

    475:     // Have an improvement for the merit function.  Calculate the function
    476:     // values for the global merit function.

    478:     printf("ip_merit: %5.4e\n", isls->merit);

    480:     Tao_ASLS_FunctionGradient(tao, x, &psi, isls->dpsi, isls);
    481:     info = isls->dpsi->Norm2(&ndpsi); CHKERRQ(info);
    482:   }

    484:   info = TaoVecDestroy(isls->dm); CHKERRQ(info);
    485:   isls->dm = TAO_NULL;

    487:   info = TaoMatDestroy(Jsub); CHKERRQ(info);
    488:   info = TaoVecDestroy(r1); CHKERRQ(info);
    489:   info = TaoVecDestroy(r2); CHKERRQ(info);
    490:   info = TaoVecDestroy(r3); CHKERRQ(info);

    492:   info = TaoIndexSetDestroy(FreeVariableSet); CHKERRQ(info);
    493:   info = TaoIndexSetDestroy(FixedVariableSet); CHKERRQ(info);
    494:   TaoFunctionReturn(0);
    495: }

    497: /* ---------------------------------------------------------- */
    501: int TaoCreate_ISILS(TAO_SOLVER tao)
    502: {
    503:   TAO_SSLS *isls;
    504:   int        info;

    506:   TaoFunctionBegin;

    508:   info = TaoNew(TAO_SSLS,&isls); CHKERRQ(info);
    509:   info = PetscLogObjectMemory(tao, sizeof(TAO_SSLS)); CHKERRQ(info);

    511:   isls->delta = 1e-10;
    512:   isls->rho = 2.1;
    513:   isls->mu_init = 0.1;

    515:   isls->identifier = 1e-5;

    517:   info=TaoSetTaoSolveRoutine(tao,TaoSolve_ISILS,(void*)isls); CHKERRQ(info);
    518:   info=TaoSetTaoSetUpDownRoutines(tao,TaoSetUp_SSLS,TaoSetDown_SSLS); CHKERRQ(info);
    519:   info=TaoSetTaoOptionsRoutine(tao,TaoSetOptions_SSLS); CHKERRQ(info);
    520:   info=TaoSetTaoViewRoutine(tao,TaoView_SSLS); CHKERRQ(info);

    522:   info = TaoCreateArmijoLineSearch(tao); CHKERRQ(info);

    524:   info = TaoSetMaximumIterates(tao,2000); CHKERRQ(info);
    525:   info = TaoSetMaximumFunctionEvaluations(tao,4000); CHKERRQ(info);

    527:   info = TaoSetTolerances(tao,0,0,0,0); CHKERRQ(info);
    528:   info = TaoSetGradientTolerances(tao,1.0e-16,0.0,0.0); CHKERRQ(info);
    529:   info = TaoSetFunctionLowerBound(tao,1.0e-8); CHKERRQ(info);

    531:   TaoFunctionReturn(0);
    532: }