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