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;
149:   TaoInt  nn, nf;
150:   TaoTerminateReason reason;
151:   TaoTruth flag;

153:   TaoFunctionBegin;

155:   // Assume that Setup has been called!

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

251:   while (1) {

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

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

261:     // Check the termination criteria

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

349:     info = FixedVariableSet->GetSize(&nf);
350:     if (nf) PetscPrintf(((PetscObject)tao)->comm,"Fixed size: %d\n", nf);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

507:   TaoFunctionBegin;

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

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

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

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

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

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

530:   TaoFunctionReturn(0);
531: }