Actual source code: bmrm.c

  1: #include <../src/tao/unconstrained/impls/bmrm/bmrm.h>

  3: static PetscErrorCode init_df_solver(TAO_DF *);
  4: static PetscErrorCode ensure_df_space(PetscInt, TAO_DF *);
  5: static PetscErrorCode destroy_df_solver(TAO_DF *);
  6: static PetscReal      phi(PetscReal *, PetscInt, PetscReal, PetscReal *, PetscReal, PetscReal *, PetscReal *, PetscReal *);
  7: static PetscInt       project(PetscInt, PetscReal *, PetscReal, PetscReal *, PetscReal *, PetscReal *, PetscReal *, PetscReal *, TAO_DF *);
  8: static PetscErrorCode solve(TAO_DF *);

 10: /*------------------------------------------------------------*/
 11: /* The main solver function

 13:    f = Remp(W)          This is what the user provides us from the application layer
 14:    So the ComputeGradient function for instance should get us back the subgradient of Remp(W)

 16:    Regularizer assumed to be L2 norm = lambda*0.5*W'W ()
 17: */

 19: static PetscErrorCode make_grad_node(Vec X, Vec_Chain **p)
 20: {
 21:   PetscNew(p);
 22:   VecDuplicate(X, &(*p)->V);
 23:   VecCopy(X, (*p)->V);
 24:   (*p)->next = NULL;
 25:   return 0;
 26: }

 28: static PetscErrorCode destroy_grad_list(Vec_Chain *head)
 29: {
 30:   Vec_Chain *p = head->next, *q;

 32:   while (p) {
 33:     q = p->next;
 34:     VecDestroy(&p->V);
 35:     PetscFree(p);
 36:     p = q;
 37:   }
 38:   head->next = NULL;
 39:   return 0;
 40: }

 42: static PetscErrorCode TaoSolve_BMRM(Tao tao)
 43: {
 44:   TAO_DF    df;
 45:   TAO_BMRM *bmrm = (TAO_BMRM *)tao->data;

 47:   /* Values and pointers to parts of the optimization problem */
 48:   PetscReal   f = 0.0;
 49:   Vec         W = tao->solution;
 50:   Vec         G = tao->gradient;
 51:   PetscReal   lambda;
 52:   PetscReal   bt;
 53:   Vec_Chain   grad_list, *tail_glist, *pgrad;
 54:   PetscInt    i;
 55:   PetscMPIInt rank;

 57:   /* Used in converged criteria check */
 58:   PetscReal reg;
 59:   PetscReal jtwt = 0.0, max_jtwt, pre_epsilon, epsilon, jw, min_jw;
 60:   PetscReal innerSolverTol;
 61:   MPI_Comm  comm;

 63:   PetscObjectGetComm((PetscObject)tao, &comm);
 64:   MPI_Comm_rank(comm, &rank);
 65:   lambda = bmrm->lambda;

 67:   /* Check Stopping Condition */
 68:   tao->step      = 1.0;
 69:   max_jtwt       = -BMRM_INFTY;
 70:   min_jw         = BMRM_INFTY;
 71:   innerSolverTol = 1.0;
 72:   epsilon        = 0.0;

 74:   if (rank == 0) {
 75:     init_df_solver(&df);
 76:     grad_list.next = NULL;
 77:     tail_glist     = &grad_list;
 78:   }

 80:   df.tol      = 1e-6;
 81:   tao->reason = TAO_CONTINUE_ITERATING;

 83:   /*-----------------Algorithm Begins------------------------*/
 84:   /* make the scatter */
 85:   VecScatterCreateToZero(W, &bmrm->scatter, &bmrm->local_w);
 86:   VecAssemblyBegin(bmrm->local_w);
 87:   VecAssemblyEnd(bmrm->local_w);

 89:   /* NOTE: In application pass the sub-gradient of Remp(W) */
 90:   TaoComputeObjectiveAndGradient(tao, W, &f, G);
 91:   TaoLogConvergenceHistory(tao, f, 1.0, 0.0, tao->ksp_its);
 92:   TaoMonitor(tao, tao->niter, f, 1.0, 0.0, tao->step);
 93:   PetscUseTypeMethod(tao, convergencetest, tao->cnvP);

 95:   while (tao->reason == TAO_CONTINUE_ITERATING) {
 96:     /* Call general purpose update function */
 97:     PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);

 99:     /* compute bt = Remp(Wt-1) - <Wt-1, At> */
100:     VecDot(W, G, &bt);
101:     bt = f - bt;

103:     /* First gather the gradient to the rank-0 node */
104:     VecScatterBegin(bmrm->scatter, G, bmrm->local_w, INSERT_VALUES, SCATTER_FORWARD);
105:     VecScatterEnd(bmrm->scatter, G, bmrm->local_w, INSERT_VALUES, SCATTER_FORWARD);

107:     /* Bring up the inner solver */
108:     if (rank == 0) {
109:       ensure_df_space(tao->niter + 1, &df);
110:       make_grad_node(bmrm->local_w, &pgrad);
111:       tail_glist->next = pgrad;
112:       tail_glist       = pgrad;

114:       df.a[tao->niter] = 1.0;
115:       df.f[tao->niter] = -bt;
116:       df.u[tao->niter] = 1.0;
117:       df.l[tao->niter] = 0.0;

119:       /* set up the Q */
120:       pgrad = grad_list.next;
121:       for (i = 0; i <= tao->niter; i++) {
123:         VecDot(pgrad->V, bmrm->local_w, &reg);
124:         df.Q[i][tao->niter] = df.Q[tao->niter][i] = reg / lambda;
125:         pgrad                                     = pgrad->next;
126:       }

128:       if (tao->niter > 0) {
129:         df.x[tao->niter] = 0.0;
130:         solve(&df);
131:       } else df.x[0] = 1.0;

133:       /* now computing Jt*(alpha_t) which should be = Jt(wt) to check convergence */
134:       jtwt = 0.0;
135:       VecSet(bmrm->local_w, 0.0);
136:       pgrad = grad_list.next;
137:       for (i = 0; i <= tao->niter; i++) {
138:         jtwt -= df.x[i] * df.f[i];
139:         VecAXPY(bmrm->local_w, -df.x[i] / lambda, pgrad->V);
140:         pgrad = pgrad->next;
141:       }

143:       VecNorm(bmrm->local_w, NORM_2, &reg);
144:       reg = 0.5 * lambda * reg * reg;
145:       jtwt -= reg;
146:     } /* end if rank == 0 */

148:     /* scatter the new W to all nodes */
149:     VecScatterBegin(bmrm->scatter, bmrm->local_w, W, INSERT_VALUES, SCATTER_REVERSE);
150:     VecScatterEnd(bmrm->scatter, bmrm->local_w, W, INSERT_VALUES, SCATTER_REVERSE);

152:     TaoComputeObjectiveAndGradient(tao, W, &f, G);

154:     MPI_Bcast(&jtwt, 1, MPIU_REAL, 0, comm);
155:     MPI_Bcast(&reg, 1, MPIU_REAL, 0, comm);

157:     jw = reg + f; /* J(w) = regularizer + Remp(w) */
158:     if (jw < min_jw) min_jw = jw;
159:     if (jtwt > max_jtwt) max_jtwt = jtwt;

161:     pre_epsilon = epsilon;
162:     epsilon     = min_jw - jtwt;

164:     if (rank == 0) {
165:       if (innerSolverTol > epsilon) innerSolverTol = epsilon;
166:       else if (innerSolverTol < 1e-7) innerSolverTol = 1e-7;

168:       /* if the annealing doesn't work well, lower the inner solver tolerance */
169:       if (pre_epsilon < epsilon) innerSolverTol *= 0.2;

171:       df.tol = innerSolverTol * 0.5;
172:     }

174:     tao->niter++;
175:     TaoLogConvergenceHistory(tao, min_jw, epsilon, 0.0, tao->ksp_its);
176:     TaoMonitor(tao, tao->niter, min_jw, epsilon, 0.0, tao->step);
177:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
178:   }

180:   /* free all the memory */
181:   if (rank == 0) {
182:     destroy_grad_list(&grad_list);
183:     destroy_df_solver(&df);
184:   }

186:   VecDestroy(&bmrm->local_w);
187:   VecScatterDestroy(&bmrm->scatter);
188:   return 0;
189: }

191: /* ---------------------------------------------------------- */

193: static PetscErrorCode TaoSetup_BMRM(Tao tao)
194: {
195:   /* Allocate some arrays */
196:   if (!tao->gradient) VecDuplicate(tao->solution, &tao->gradient);
197:   return 0;
198: }

200: /*------------------------------------------------------------*/
201: static PetscErrorCode TaoDestroy_BMRM(Tao tao)
202: {
203:   PetscFree(tao->data);
204:   return 0;
205: }

207: static PetscErrorCode TaoSetFromOptions_BMRM(Tao tao, PetscOptionItems *PetscOptionsObject)
208: {
209:   TAO_BMRM *bmrm = (TAO_BMRM *)tao->data;

211:   PetscOptionsHeadBegin(PetscOptionsObject, "BMRM for regularized risk minimization");
212:   PetscOptionsReal("-tao_bmrm_lambda", "regulariser weight", "", 100, &bmrm->lambda, NULL);
213:   PetscOptionsHeadEnd();
214:   return 0;
215: }

217: /*------------------------------------------------------------*/
218: static PetscErrorCode TaoView_BMRM(Tao tao, PetscViewer viewer)
219: {
220:   PetscBool isascii;

222:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
223:   if (isascii) {
224:     PetscViewerASCIIPushTab(viewer);
225:     PetscViewerASCIIPopTab(viewer);
226:   }
227:   return 0;
228: }

230: /*------------------------------------------------------------*/
231: /*MC
232:   TAOBMRM - bundle method for regularized risk minimization

234:   Options Database Keys:
235: . - tao_bmrm_lambda - regulariser weight

237:   Level: beginner
238: M*/

240: PETSC_EXTERN PetscErrorCode TaoCreate_BMRM(Tao tao)
241: {
242:   TAO_BMRM *bmrm;

244:   tao->ops->setup          = TaoSetup_BMRM;
245:   tao->ops->solve          = TaoSolve_BMRM;
246:   tao->ops->view           = TaoView_BMRM;
247:   tao->ops->setfromoptions = TaoSetFromOptions_BMRM;
248:   tao->ops->destroy        = TaoDestroy_BMRM;

250:   PetscNew(&bmrm);
251:   bmrm->lambda = 1.0;
252:   tao->data    = (void *)bmrm;

254:   /* Override default settings (unless already changed) */
255:   if (!tao->max_it_changed) tao->max_it = 2000;
256:   if (!tao->max_funcs_changed) tao->max_funcs = 4000;
257:   if (!tao->gatol_changed) tao->gatol = 1.0e-12;
258:   if (!tao->grtol_changed) tao->grtol = 1.0e-12;

260:   return 0;
261: }

263: PetscErrorCode init_df_solver(TAO_DF *df)
264: {
265:   PetscInt i, n = INCRE_DIM;

267:   /* default values */
268:   df->maxProjIter = 200;
269:   df->maxPGMIter  = 300000;
270:   df->b           = 1.0;

272:   /* memory space required by Dai-Fletcher */
273:   df->cur_num_cp = n;
274:   PetscMalloc1(n, &df->f);
275:   PetscMalloc1(n, &df->a);
276:   PetscMalloc1(n, &df->l);
277:   PetscMalloc1(n, &df->u);
278:   PetscMalloc1(n, &df->x);
279:   PetscMalloc1(n, &df->Q);

281:   for (i = 0; i < n; i++) PetscMalloc1(n, &df->Q[i]);

283:   PetscMalloc1(n, &df->g);
284:   PetscMalloc1(n, &df->y);
285:   PetscMalloc1(n, &df->tempv);
286:   PetscMalloc1(n, &df->d);
287:   PetscMalloc1(n, &df->Qd);
288:   PetscMalloc1(n, &df->t);
289:   PetscMalloc1(n, &df->xplus);
290:   PetscMalloc1(n, &df->tplus);
291:   PetscMalloc1(n, &df->sk);
292:   PetscMalloc1(n, &df->yk);

294:   PetscMalloc1(n, &df->ipt);
295:   PetscMalloc1(n, &df->ipt2);
296:   PetscMalloc1(n, &df->uv);
297:   return 0;
298: }

300: PetscErrorCode ensure_df_space(PetscInt dim, TAO_DF *df)
301: {
302:   PetscReal *tmp, **tmp_Q;
303:   PetscInt   i, n, old_n;

305:   df->dim = dim;
306:   if (dim <= df->cur_num_cp) return 0;

308:   old_n = df->cur_num_cp;
309:   df->cur_num_cp += INCRE_DIM;
310:   n = df->cur_num_cp;

312:   /* memory space required by dai-fletcher */
313:   PetscMalloc1(n, &tmp);
314:   PetscArraycpy(tmp, df->f, old_n);
315:   PetscFree(df->f);
316:   df->f = tmp;

318:   PetscMalloc1(n, &tmp);
319:   PetscArraycpy(tmp, df->a, old_n);
320:   PetscFree(df->a);
321:   df->a = tmp;

323:   PetscMalloc1(n, &tmp);
324:   PetscArraycpy(tmp, df->l, old_n);
325:   PetscFree(df->l);
326:   df->l = tmp;

328:   PetscMalloc1(n, &tmp);
329:   PetscArraycpy(tmp, df->u, old_n);
330:   PetscFree(df->u);
331:   df->u = tmp;

333:   PetscMalloc1(n, &tmp);
334:   PetscArraycpy(tmp, df->x, old_n);
335:   PetscFree(df->x);
336:   df->x = tmp;

338:   PetscMalloc1(n, &tmp_Q);
339:   for (i = 0; i < n; i++) {
340:     PetscMalloc1(n, &tmp_Q[i]);
341:     if (i < old_n) {
342:       PetscArraycpy(tmp_Q[i], df->Q[i], old_n);
343:       PetscFree(df->Q[i]);
344:     }
345:   }

347:   PetscFree(df->Q);
348:   df->Q = tmp_Q;

350:   PetscFree(df->g);
351:   PetscMalloc1(n, &df->g);

353:   PetscFree(df->y);
354:   PetscMalloc1(n, &df->y);

356:   PetscFree(df->tempv);
357:   PetscMalloc1(n, &df->tempv);

359:   PetscFree(df->d);
360:   PetscMalloc1(n, &df->d);

362:   PetscFree(df->Qd);
363:   PetscMalloc1(n, &df->Qd);

365:   PetscFree(df->t);
366:   PetscMalloc1(n, &df->t);

368:   PetscFree(df->xplus);
369:   PetscMalloc1(n, &df->xplus);

371:   PetscFree(df->tplus);
372:   PetscMalloc1(n, &df->tplus);

374:   PetscFree(df->sk);
375:   PetscMalloc1(n, &df->sk);

377:   PetscFree(df->yk);
378:   PetscMalloc1(n, &df->yk);

380:   PetscFree(df->ipt);
381:   PetscMalloc1(n, &df->ipt);

383:   PetscFree(df->ipt2);
384:   PetscMalloc1(n, &df->ipt2);

386:   PetscFree(df->uv);
387:   PetscMalloc1(n, &df->uv);
388:   return 0;
389: }

391: PetscErrorCode destroy_df_solver(TAO_DF *df)
392: {
393:   PetscInt i;

395:   PetscFree(df->f);
396:   PetscFree(df->a);
397:   PetscFree(df->l);
398:   PetscFree(df->u);
399:   PetscFree(df->x);

401:   for (i = 0; i < df->cur_num_cp; i++) PetscFree(df->Q[i]);
402:   PetscFree(df->Q);
403:   PetscFree(df->ipt);
404:   PetscFree(df->ipt2);
405:   PetscFree(df->uv);
406:   PetscFree(df->g);
407:   PetscFree(df->y);
408:   PetscFree(df->tempv);
409:   PetscFree(df->d);
410:   PetscFree(df->Qd);
411:   PetscFree(df->t);
412:   PetscFree(df->xplus);
413:   PetscFree(df->tplus);
414:   PetscFree(df->sk);
415:   PetscFree(df->yk);
416:   return 0;
417: }

419: /* Piecewise linear monotone target function for the Dai-Fletcher projector */
420: PetscReal phi(PetscReal *x, PetscInt n, PetscReal lambda, PetscReal *a, PetscReal b, PetscReal *c, PetscReal *l, PetscReal *u)
421: {
422:   PetscReal r = 0.0;
423:   PetscInt  i;

425:   for (i = 0; i < n; i++) {
426:     x[i] = -c[i] + lambda * a[i];
427:     if (x[i] > u[i]) x[i] = u[i];
428:     else if (x[i] < l[i]) x[i] = l[i];
429:     r += a[i] * x[i];
430:   }
431:   return r - b;
432: }

434: /** Modified Dai-Fletcher QP projector solves the problem:
435:  *
436:  *      minimise  0.5*x'*x - c'*x
437:  *      subj to   a'*x = b
438:  *                l \leq x \leq u
439:  *
440:  *  \param c The point to be projected onto feasible set
441:  */
442: PetscInt project(PetscInt n, PetscReal *a, PetscReal b, PetscReal *c, PetscReal *l, PetscReal *u, PetscReal *x, PetscReal *lam_ext, TAO_DF *df)
443: {
444:   PetscReal lambda, lambdal, lambdau, dlambda, lambda_new;
445:   PetscReal r, rl, ru, s;
446:   PetscInt  innerIter;
447:   PetscBool nonNegativeSlack = PETSC_FALSE;

449:   *lam_ext  = 0;
450:   lambda    = 0;
451:   dlambda   = 0.5;
452:   innerIter = 1;

454:   /*  \phi(x;lambda) := 0.5*x'*x + c'*x - lambda*(a'*x-b)
455:    *
456:    *  Optimality conditions for \phi:
457:    *
458:    *  1. lambda   <= 0
459:    *  2. r        <= 0
460:    *  3. r*lambda == 0
461:    */

463:   /* Bracketing Phase */
464:   r = phi(x, n, lambda, a, b, c, l, u);

466:   if (nonNegativeSlack) {
467:     /* inequality constraint, i.e., with \xi >= 0 constraint */
468:     if (r < TOL_R) return 0;
469:   } else {
470:     /* equality constraint ,i.e., without \xi >= 0 constraint */
471:     if (PetscAbsReal(r) < TOL_R) return 0;
472:   }

474:   if (r < 0.0) {
475:     lambdal = lambda;
476:     rl      = r;
477:     lambda  = lambda + dlambda;
478:     r       = phi(x, n, lambda, a, b, c, l, u);
479:     while (r < 0.0 && dlambda < BMRM_INFTY) {
480:       lambdal = lambda;
481:       s       = rl / r - 1.0;
482:       if (s < 0.1) s = 0.1;
483:       dlambda = dlambda + dlambda / s;
484:       lambda  = lambda + dlambda;
485:       rl      = r;
486:       r       = phi(x, n, lambda, a, b, c, l, u);
487:     }
488:     lambdau = lambda;
489:     ru      = r;
490:   } else {
491:     lambdau = lambda;
492:     ru      = r;
493:     lambda  = lambda - dlambda;
494:     r       = phi(x, n, lambda, a, b, c, l, u);
495:     while (r > 0.0 && dlambda > -BMRM_INFTY) {
496:       lambdau = lambda;
497:       s       = ru / r - 1.0;
498:       if (s < 0.1) s = 0.1;
499:       dlambda = dlambda + dlambda / s;
500:       lambda  = lambda - dlambda;
501:       ru      = r;
502:       r       = phi(x, n, lambda, a, b, c, l, u);
503:     }
504:     lambdal = lambda;
505:     rl      = r;
506:   }


510:   if (ru == 0) return innerIter;

512:   /* Secant Phase */
513:   s       = 1.0 - rl / ru;
514:   dlambda = dlambda / s;
515:   lambda  = lambdau - dlambda;
516:   r       = phi(x, n, lambda, a, b, c, l, u);

518:   while (PetscAbsReal(r) > TOL_R && dlambda > TOL_LAM * (1.0 + PetscAbsReal(lambda)) && innerIter < df->maxProjIter) {
519:     innerIter++;
520:     if (r > 0.0) {
521:       if (s <= 2.0) {
522:         lambdau = lambda;
523:         ru      = r;
524:         s       = 1.0 - rl / ru;
525:         dlambda = (lambdau - lambdal) / s;
526:         lambda  = lambdau - dlambda;
527:       } else {
528:         s = ru / r - 1.0;
529:         if (s < 0.1) s = 0.1;
530:         dlambda    = (lambdau - lambda) / s;
531:         lambda_new = 0.75 * lambdal + 0.25 * lambda;
532:         if (lambda_new < (lambda - dlambda)) lambda_new = lambda - dlambda;
533:         lambdau = lambda;
534:         ru      = r;
535:         lambda  = lambda_new;
536:         s       = (lambdau - lambdal) / (lambdau - lambda);
537:       }
538:     } else {
539:       if (s >= 2.0) {
540:         lambdal = lambda;
541:         rl      = r;
542:         s       = 1.0 - rl / ru;
543:         dlambda = (lambdau - lambdal) / s;
544:         lambda  = lambdau - dlambda;
545:       } else {
546:         s = rl / r - 1.0;
547:         if (s < 0.1) s = 0.1;
548:         dlambda    = (lambda - lambdal) / s;
549:         lambda_new = 0.75 * lambdau + 0.25 * lambda;
550:         if (lambda_new > (lambda + dlambda)) lambda_new = lambda + dlambda;
551:         lambdal = lambda;
552:         rl      = r;
553:         lambda  = lambda_new;
554:         s       = (lambdau - lambdal) / (lambdau - lambda);
555:       }
556:     }
557:     r = phi(x, n, lambda, a, b, c, l, u);
558:   }

560:   *lam_ext = lambda;
561:   if (innerIter >= df->maxProjIter) PetscInfo(NULL, "WARNING: DaiFletcher max iterations\n");
562:   return innerIter;
563: }

565: PetscErrorCode solve(TAO_DF *df)
566: {
567:   PetscInt    i, j, innerIter, it, it2, luv, info;
568:   PetscReal   gd, max, ak, bk, akold, bkold, lamnew, alpha, kktlam = 0.0, lam_ext;
569:   PetscReal   DELTAsv, ProdDELTAsv;
570:   PetscReal   c, *tempQ;
571:   PetscReal  *x = df->x, *a = df->a, b = df->b, *l = df->l, *u = df->u, tol = df->tol;
572:   PetscReal  *tempv = df->tempv, *y = df->y, *g = df->g, *d = df->d, *Qd = df->Qd;
573:   PetscReal  *xplus = df->xplus, *tplus = df->tplus, *sk = df->sk, *yk = df->yk;
574:   PetscReal **Q = df->Q, *f = df->f, *t = df->t;
575:   PetscInt    dim = df->dim, *ipt = df->ipt, *ipt2 = df->ipt2, *uv = df->uv;

577:   /* variables for the adaptive nonmonotone linesearch */
578:   PetscInt  L, llast;
579:   PetscReal fr, fbest, fv, fc, fv0;

581:   c = BMRM_INFTY;

583:   DELTAsv = EPS_SV;
584:   if (tol <= 1.0e-5 || dim <= 20) ProdDELTAsv = 0.0F;
585:   else ProdDELTAsv = EPS_SV;

587:   for (i = 0; i < dim; i++) tempv[i] = -x[i];

589:   lam_ext = 0.0;

591:   /* Project the initial solution */
592:   project(dim, a, b, tempv, l, u, x, &lam_ext, df);

594:   /* Compute gradient
595:      g = Q*x + f; */

597:   it = 0;
598:   for (i = 0; i < dim; i++) {
599:     if (PetscAbsReal(x[i]) > ProdDELTAsv) ipt[it++] = i;
600:   }

602:   PetscArrayzero(t, dim);
603:   for (i = 0; i < it; i++) {
604:     tempQ = Q[ipt[i]];
605:     for (j = 0; j < dim; j++) t[j] += (tempQ[j] * x[ipt[i]]);
606:   }
607:   for (i = 0; i < dim; i++) g[i] = t[i] + f[i];

609:   /* y = -(x_{k} - g_{k}) */
610:   for (i = 0; i < dim; i++) y[i] = g[i] - x[i];

612:   /* Project x_{k} - g_{k} */
613:   project(dim, a, b, y, l, u, tempv, &lam_ext, df);

615:   /* y = P(x_{k} - g_{k}) - x_{k} */
616:   max = ALPHA_MIN;
617:   for (i = 0; i < dim; i++) {
618:     y[i] = tempv[i] - x[i];
619:     if (PetscAbsReal(y[i]) > max) max = PetscAbsReal(y[i]);
620:   }

622:   if (max < tol * 1e-3) return 0;

624:   alpha = 1.0 / max;

626:   /* fv0 = f(x_{0}). Recall t = Q x_{k}  */
627:   fv0 = 0.0;
628:   for (i = 0; i < dim; i++) fv0 += x[i] * (0.5 * t[i] + f[i]);

630:   /* adaptive nonmonotone linesearch */
631:   L     = 2;
632:   fr    = ALPHA_MAX;
633:   fbest = fv0;
634:   fc    = fv0;
635:   llast = 0;
636:   akold = bkold = 0.0;

638:   /*     Iterator begins     */
639:   for (innerIter = 1; innerIter <= df->maxPGMIter; innerIter++) {
640:     /* tempv = -(x_{k} - alpha*g_{k}) */
641:     for (i = 0; i < dim; i++) tempv[i] = alpha * g[i] - x[i];

643:     /* Project x_{k} - alpha*g_{k} */
644:     project(dim, a, b, tempv, l, u, y, &lam_ext, df);

646:     /* gd = \inner{d_{k}}{g_{k}}
647:         d = P(x_{k} - alpha*g_{k}) - x_{k}
648:     */
649:     gd = 0.0;
650:     for (i = 0; i < dim; i++) {
651:       d[i] = y[i] - x[i];
652:       gd += d[i] * g[i];
653:     }

655:     /* Gradient computation  */

657:     /* compute Qd = Q*d  or  Qd = Q*y - t depending on their sparsity */

659:     it = it2 = 0;
660:     for (i = 0; i < dim; i++) {
661:       if (PetscAbsReal(d[i]) > (ProdDELTAsv * 1.0e-2)) ipt[it++] = i;
662:     }
663:     for (i = 0; i < dim; i++) {
664:       if (PetscAbsReal(y[i]) > ProdDELTAsv) ipt2[it2++] = i;
665:     }

667:     PetscArrayzero(Qd, dim);
668:     /* compute Qd = Q*d */
669:     if (it < it2) {
670:       for (i = 0; i < it; i++) {
671:         tempQ = Q[ipt[i]];
672:         for (j = 0; j < dim; j++) Qd[j] += (tempQ[j] * d[ipt[i]]);
673:       }
674:     } else { /* compute Qd = Q*y-t */
675:       for (i = 0; i < it2; i++) {
676:         tempQ = Q[ipt2[i]];
677:         for (j = 0; j < dim; j++) Qd[j] += (tempQ[j] * y[ipt2[i]]);
678:       }
679:       for (j = 0; j < dim; j++) Qd[j] -= t[j];
680:     }

682:     /* ak = inner{d_{k}}{d_{k}} */
683:     ak = 0.0;
684:     for (i = 0; i < dim; i++) ak += d[i] * d[i];

686:     bk = 0.0;
687:     for (i = 0; i < dim; i++) bk += d[i] * Qd[i];

689:     if (bk > EPS * ak && gd < 0.0) lamnew = -gd / bk;
690:     else lamnew = 1.0;

692:     /* fv is computing f(x_{k} + d_{k}) */
693:     fv = 0.0;
694:     for (i = 0; i < dim; i++) {
695:       xplus[i] = x[i] + d[i];
696:       tplus[i] = t[i] + Qd[i];
697:       fv += xplus[i] * (0.5 * tplus[i] + f[i]);
698:     }

700:     /* fr is fref */
701:     if ((innerIter == 1 && fv >= fv0) || (innerIter > 1 && fv >= fr)) {
702:       fv = 0.0;
703:       for (i = 0; i < dim; i++) {
704:         xplus[i] = x[i] + lamnew * d[i];
705:         tplus[i] = t[i] + lamnew * Qd[i];
706:         fv += xplus[i] * (0.5 * tplus[i] + f[i]);
707:       }
708:     }

710:     for (i = 0; i < dim; i++) {
711:       sk[i] = xplus[i] - x[i];
712:       yk[i] = tplus[i] - t[i];
713:       x[i]  = xplus[i];
714:       t[i]  = tplus[i];
715:       g[i]  = t[i] + f[i];
716:     }

718:     /* update the line search control parameters */
719:     if (fv < fbest) {
720:       fbest = fv;
721:       fc    = fv;
722:       llast = 0;
723:     } else {
724:       fc = (fc > fv ? fc : fv);
725:       llast++;
726:       if (llast == L) {
727:         fr    = fc;
728:         fc    = fv;
729:         llast = 0;
730:       }
731:     }

733:     ak = bk = 0.0;
734:     for (i = 0; i < dim; i++) {
735:       ak += sk[i] * sk[i];
736:       bk += sk[i] * yk[i];
737:     }

739:     if (bk <= EPS * ak) alpha = ALPHA_MAX;
740:     else {
741:       if (bkold < EPS * akold) alpha = ak / bk;
742:       else alpha = (akold + ak) / (bkold + bk);

744:       if (alpha > ALPHA_MAX) alpha = ALPHA_MAX;
745:       else if (alpha < ALPHA_MIN) alpha = ALPHA_MIN;
746:     }

748:     akold = ak;
749:     bkold = bk;

751:     /* stopping criterion based on KKT conditions */
752:     /* at optimal, gradient of lagrangian w.r.t. x is zero */

754:     bk = 0.0;
755:     for (i = 0; i < dim; i++) bk += x[i] * x[i];

757:     if (PetscSqrtReal(ak) < tol * 10 * PetscSqrtReal(bk)) {
758:       it     = 0;
759:       luv    = 0;
760:       kktlam = 0.0;
761:       for (i = 0; i < dim; i++) {
762:         /* x[i] is active hence lagrange multipliers for box constraints
763:                 are zero. The lagrange multiplier for ineq. const. is then
764:                 defined as below
765:         */
766:         if ((x[i] > DELTAsv) && (x[i] < c - DELTAsv)) {
767:           ipt[it++] = i;
768:           kktlam    = kktlam - a[i] * g[i];
769:         } else uv[luv++] = i;
770:       }

772:       if (it == 0 && PetscSqrtReal(ak) < tol * 0.5 * PetscSqrtReal(bk)) return 0;
773:       else {
774:         kktlam = kktlam / it;
775:         info   = 1;
776:         for (i = 0; i < it; i++) {
777:           if (PetscAbsReal(a[ipt[i]] * g[ipt[i]] + kktlam) > tol) {
778:             info = 0;
779:             break;
780:           }
781:         }
782:         if (info == 1) {
783:           for (i = 0; i < luv; i++) {
784:             if (x[uv[i]] <= DELTAsv) {
785:               /* x[i] == lower bound, hence, lagrange multiplier (say, beta) for lower bound may
786:                      not be zero. So, the gradient without beta is > 0
787:               */
788:               if (g[uv[i]] + kktlam * a[uv[i]] < -tol) {
789:                 info = 0;
790:                 break;
791:               }
792:             } else {
793:               /* x[i] == upper bound, hence, lagrange multiplier (say, eta) for upper bound may
794:                      not be zero. So, the gradient without eta is < 0
795:               */
796:               if (g[uv[i]] + kktlam * a[uv[i]] > tol) {
797:                 info = 0;
798:                 break;
799:               }
800:             }
801:           }
802:         }

804:         if (info == 1) return 0;
805:       }
806:     }
807:   }
808:   return 0;
809: }