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 *);

  9: static PetscErrorCode solve(TAO_DF *df)
 10: {
 11:   PetscInt    i, j, innerIter, it, it2, luv, info;
 12:   PetscReal   gd, max, ak, bk, akold, bkold, lamnew, alpha, kktlam = 0.0, lam_ext;
 13:   PetscReal   DELTAsv, ProdDELTAsv;
 14:   PetscReal   c, *tempQ;
 15:   PetscReal  *x = df->x, *a = df->a, b = df->b, *l = df->l, *u = df->u, tol = df->tol;
 16:   PetscReal  *tempv = df->tempv, *y = df->y, *g = df->g, *d = df->d, *Qd = df->Qd;
 17:   PetscReal  *xplus = df->xplus, *tplus = df->tplus, *sk = df->sk, *yk = df->yk;
 18:   PetscReal **Q = df->Q, *f = df->f, *t = df->t;
 19:   PetscInt    dim = df->dim, *ipt = df->ipt, *ipt2 = df->ipt2, *uv = df->uv;

 21:   /* variables for the adaptive nonmonotone linesearch */
 22:   PetscInt  L, llast;
 23:   PetscReal fr, fbest, fv, fc, fv0;

 25:   c = BMRM_INFTY;

 27:   DELTAsv = EPS_SV;
 28:   if (tol <= 1.0e-5 || dim <= 20) ProdDELTAsv = 0.0F;
 29:   else ProdDELTAsv = EPS_SV;

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

 33:   lam_ext = 0.0;

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

 38:   /* Compute gradient
 39:      g = Q*x + f; */

 41:   it = 0;
 42:   for (i = 0; i < dim; i++) {
 43:     if (PetscAbsReal(x[i]) > ProdDELTAsv) ipt[it++] = i;
 44:   }

 46:   PetscCall(PetscArrayzero(t, dim));
 47:   for (i = 0; i < it; i++) {
 48:     tempQ = Q[ipt[i]];
 49:     for (j = 0; j < dim; j++) t[j] += (tempQ[j] * x[ipt[i]]);
 50:   }
 51:   for (i = 0; i < dim; i++) g[i] = t[i] + f[i];

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

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

 59:   /* y = P(x_{k} - g_{k}) - x_{k} */
 60:   max = ALPHA_MIN;
 61:   for (i = 0; i < dim; i++) {
 62:     y[i] = tempv[i] - x[i];
 63:     if (PetscAbsReal(y[i]) > max) max = PetscAbsReal(y[i]);
 64:   }

 66:   if (max < tol * 1e-3) return PETSC_SUCCESS;

 68:   alpha = 1.0 / max;

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

 74:   /* adaptive nonmonotone linesearch */
 75:   L     = 2;
 76:   fr    = ALPHA_MAX;
 77:   fbest = fv0;
 78:   fc    = fv0;
 79:   llast = 0;
 80:   akold = bkold = 0.0;

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

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

 90:     /* gd = \inner{d_{k}}{g_{k}}
 91:         d = P(x_{k} - alpha*g_{k}) - x_{k}
 92:     */
 93:     gd = 0.0;
 94:     for (i = 0; i < dim; i++) {
 95:       d[i] = y[i] - x[i];
 96:       gd += d[i] * g[i];
 97:     }

 99:     /* Gradient computation  */

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

103:     it = it2 = 0;
104:     for (i = 0; i < dim; i++) {
105:       if (PetscAbsReal(d[i]) > (ProdDELTAsv * 1.0e-2)) ipt[it++] = i;
106:     }
107:     for (i = 0; i < dim; i++) {
108:       if (PetscAbsReal(y[i]) > ProdDELTAsv) ipt2[it2++] = i;
109:     }

111:     PetscCall(PetscArrayzero(Qd, dim));
112:     /* compute Qd = Q*d */
113:     if (it < it2) {
114:       for (i = 0; i < it; i++) {
115:         tempQ = Q[ipt[i]];
116:         for (j = 0; j < dim; j++) Qd[j] += (tempQ[j] * d[ipt[i]]);
117:       }
118:     } else { /* compute Qd = Q*y-t */
119:       for (i = 0; i < it2; i++) {
120:         tempQ = Q[ipt2[i]];
121:         for (j = 0; j < dim; j++) Qd[j] += (tempQ[j] * y[ipt2[i]]);
122:       }
123:       for (j = 0; j < dim; j++) Qd[j] -= t[j];
124:     }

126:     /* ak = inner{d_{k}}{d_{k}} */
127:     ak = 0.0;
128:     for (i = 0; i < dim; i++) ak += d[i] * d[i];

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

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

136:     /* fv is computing f(x_{k} + d_{k}) */
137:     fv = 0.0;
138:     for (i = 0; i < dim; i++) {
139:       xplus[i] = x[i] + d[i];
140:       tplus[i] = t[i] + Qd[i];
141:       fv += xplus[i] * (0.5 * tplus[i] + f[i]);
142:     }

144:     /* fr is fref */
145:     if ((innerIter == 1 && fv >= fv0) || (innerIter > 1 && fv >= fr)) {
146:       fv = 0.0;
147:       for (i = 0; i < dim; i++) {
148:         xplus[i] = x[i] + lamnew * d[i];
149:         tplus[i] = t[i] + lamnew * Qd[i];
150:         fv += xplus[i] * (0.5 * tplus[i] + f[i]);
151:       }
152:     }

154:     for (i = 0; i < dim; i++) {
155:       sk[i] = xplus[i] - x[i];
156:       yk[i] = tplus[i] - t[i];
157:       x[i]  = xplus[i];
158:       t[i]  = tplus[i];
159:       g[i]  = t[i] + f[i];
160:     }

162:     /* update the line search control parameters */
163:     if (fv < fbest) {
164:       fbest = fv;
165:       fc    = fv;
166:       llast = 0;
167:     } else {
168:       fc = (fc > fv ? fc : fv);
169:       llast++;
170:       if (llast == L) {
171:         fr    = fc;
172:         fc    = fv;
173:         llast = 0;
174:       }
175:     }

177:     ak = bk = 0.0;
178:     for (i = 0; i < dim; i++) {
179:       ak += sk[i] * sk[i];
180:       bk += sk[i] * yk[i];
181:     }

183:     if (bk <= EPS * ak) alpha = ALPHA_MAX;
184:     else {
185:       if (bkold < EPS * akold) alpha = ak / bk;
186:       else alpha = (akold + ak) / (bkold + bk);

188:       if (alpha > ALPHA_MAX) alpha = ALPHA_MAX;
189:       else if (alpha < ALPHA_MIN) alpha = ALPHA_MIN;
190:     }

192:     akold = ak;
193:     bkold = bk;

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

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

201:     if (PetscSqrtReal(ak) < tol * 10 * PetscSqrtReal(bk)) {
202:       it     = 0;
203:       luv    = 0;
204:       kktlam = 0.0;
205:       for (i = 0; i < dim; i++) {
206:         /* x[i] is active hence lagrange multipliers for box constraints
207:                 are zero. The lagrange multiplier for ineq. const. is then
208:                 defined as below
209:         */
210:         if ((x[i] > DELTAsv) && (x[i] < c - DELTAsv)) {
211:           ipt[it++] = i;
212:           kktlam    = kktlam - a[i] * g[i];
213:         } else uv[luv++] = i;
214:       }

216:       if (it == 0 && PetscSqrtReal(ak) < tol * 0.5 * PetscSqrtReal(bk)) return PETSC_SUCCESS;
217:       else {
218:         kktlam = kktlam / it;
219:         info   = 1;
220:         for (i = 0; i < it; i++) {
221:           if (PetscAbsReal(a[ipt[i]] * g[ipt[i]] + kktlam) > tol) {
222:             info = 0;
223:             break;
224:           }
225:         }
226:         if (info == 1) {
227:           for (i = 0; i < luv; i++) {
228:             if (x[uv[i]] <= DELTAsv) {
229:               /* x[i] == lower bound, hence, lagrange multiplier (say, beta) for lower bound may
230:                      not be zero. So, the gradient without beta is > 0
231:               */
232:               if (g[uv[i]] + kktlam * a[uv[i]] < -tol) {
233:                 info = 0;
234:                 break;
235:               }
236:             } else {
237:               /* x[i] == upper bound, hence, lagrange multiplier (say, eta) for upper bound may
238:                      not be zero. So, the gradient without eta is < 0
239:               */
240:               if (g[uv[i]] + kktlam * a[uv[i]] > tol) {
241:                 info = 0;
242:                 break;
243:               }
244:             }
245:           }
246:         }

248:         if (info == 1) return PETSC_SUCCESS;
249:       }
250:     }
251:   }
252:   return PETSC_SUCCESS;
253: }

255: /* The main solver function

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

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

263: static PetscErrorCode make_grad_node(Vec X, Vec_Chain **p)
264: {
265:   PetscFunctionBegin;
266:   PetscCall(PetscNew(p));
267:   PetscCall(VecDuplicate(X, &(*p)->V));
268:   PetscCall(VecCopy(X, (*p)->V));
269:   (*p)->next = NULL;
270:   PetscFunctionReturn(PETSC_SUCCESS);
271: }

273: static PetscErrorCode destroy_grad_list(Vec_Chain *head)
274: {
275:   Vec_Chain *p = head->next, *q;

277:   PetscFunctionBegin;
278:   while (p) {
279:     q = p->next;
280:     PetscCall(VecDestroy(&p->V));
281:     PetscCall(PetscFree(p));
282:     p = q;
283:   }
284:   head->next = NULL;
285:   PetscFunctionReturn(PETSC_SUCCESS);
286: }

288: static PetscErrorCode TaoSolve_BMRM(Tao tao)
289: {
290:   TAO_DF    df;
291:   TAO_BMRM *bmrm = (TAO_BMRM *)tao->data;

293:   /* Values and pointers to parts of the optimization problem */
294:   PetscReal   f = 0.0;
295:   Vec         W = tao->solution;
296:   Vec         G = tao->gradient;
297:   PetscReal   lambda;
298:   PetscReal   bt;
299:   Vec_Chain   grad_list, *tail_glist, *pgrad;
300:   PetscInt    i;
301:   PetscMPIInt rank;

303:   /* Used in converged criteria check */
304:   PetscReal reg;
305:   PetscReal jtwt = 0.0, max_jtwt, pre_epsilon, epsilon, jw, min_jw;
306:   PetscReal innerSolverTol;
307:   MPI_Comm  comm;

309:   PetscFunctionBegin;
310:   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
311:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
312:   lambda = bmrm->lambda;

314:   /* Check Stopping Condition */
315:   tao->step      = 1.0;
316:   max_jtwt       = -BMRM_INFTY;
317:   min_jw         = BMRM_INFTY;
318:   innerSolverTol = 1.0;
319:   epsilon        = 0.0;

321:   if (rank == 0) {
322:     PetscCall(init_df_solver(&df));
323:     grad_list.next = NULL;
324:     tail_glist     = &grad_list;
325:   }

327:   df.tol      = 1e-6;
328:   tao->reason = TAO_CONTINUE_ITERATING;

330:   /*-----------------Algorithm Begins------------------------*/
331:   /* make the scatter */
332:   PetscCall(VecScatterCreateToZero(W, &bmrm->scatter, &bmrm->local_w));
333:   PetscCall(VecAssemblyBegin(bmrm->local_w));
334:   PetscCall(VecAssemblyEnd(bmrm->local_w));

336:   /* NOTE: In application pass the sub-gradient of Remp(W) */
337:   PetscCall(TaoComputeObjectiveAndGradient(tao, W, &f, G));
338:   PetscCall(TaoLogConvergenceHistory(tao, f, 1.0, 0.0, tao->ksp_its));
339:   PetscCall(TaoMonitor(tao, tao->niter, f, 1.0, 0.0, tao->step));
340:   PetscUseTypeMethod(tao, convergencetest, tao->cnvP);

342:   while (tao->reason == TAO_CONTINUE_ITERATING) {
343:     /* Call general purpose update function */
344:     PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);

346:     /* compute bt = Remp(Wt-1) - <Wt-1, At> */
347:     PetscCall(VecDot(W, G, &bt));
348:     bt = f - bt;

350:     /* First gather the gradient to the rank-0 node */
351:     PetscCall(VecScatterBegin(bmrm->scatter, G, bmrm->local_w, INSERT_VALUES, SCATTER_FORWARD));
352:     PetscCall(VecScatterEnd(bmrm->scatter, G, bmrm->local_w, INSERT_VALUES, SCATTER_FORWARD));

354:     /* Bring up the inner solver */
355:     if (rank == 0) {
356:       PetscCall(ensure_df_space(tao->niter + 1, &df));
357:       PetscCall(make_grad_node(bmrm->local_w, &pgrad));
358:       tail_glist->next = pgrad;
359:       tail_glist       = pgrad;

361:       df.a[tao->niter] = 1.0;
362:       df.f[tao->niter] = -bt;
363:       df.u[tao->niter] = 1.0;
364:       df.l[tao->niter] = 0.0;

366:       /* set up the Q */
367:       pgrad = grad_list.next;
368:       for (i = 0; i <= tao->niter; i++) {
369:         PetscCheck(pgrad, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Assert that there are at least tao->niter+1 pgrad available");
370:         PetscCall(VecDot(pgrad->V, bmrm->local_w, &reg));
371:         df.Q[i][tao->niter] = df.Q[tao->niter][i] = reg / lambda;
372:         pgrad                                     = pgrad->next;
373:       }

375:       if (tao->niter > 0) {
376:         df.x[tao->niter] = 0.0;
377:         PetscCall(solve(&df));
378:       } else df.x[0] = 1.0;

380:       /* now computing Jt*(alpha_t) which should be = Jt(wt) to check convergence */
381:       jtwt = 0.0;
382:       PetscCall(VecSet(bmrm->local_w, 0.0));
383:       pgrad = grad_list.next;
384:       for (i = 0; i <= tao->niter; i++) {
385:         jtwt -= df.x[i] * df.f[i];
386:         PetscCall(VecAXPY(bmrm->local_w, -df.x[i] / lambda, pgrad->V));
387:         pgrad = pgrad->next;
388:       }

390:       PetscCall(VecNorm(bmrm->local_w, NORM_2, &reg));
391:       reg = 0.5 * lambda * reg * reg;
392:       jtwt -= reg;
393:     } /* end if rank == 0 */

395:     /* scatter the new W to all nodes */
396:     PetscCall(VecScatterBegin(bmrm->scatter, bmrm->local_w, W, INSERT_VALUES, SCATTER_REVERSE));
397:     PetscCall(VecScatterEnd(bmrm->scatter, bmrm->local_w, W, INSERT_VALUES, SCATTER_REVERSE));

399:     PetscCall(TaoComputeObjectiveAndGradient(tao, W, &f, G));

401:     PetscCallMPI(MPI_Bcast(&jtwt, 1, MPIU_REAL, 0, comm));
402:     PetscCallMPI(MPI_Bcast(&reg, 1, MPIU_REAL, 0, comm));

404:     jw = reg + f; /* J(w) = regularizer + Remp(w) */
405:     if (jw < min_jw) min_jw = jw;
406:     if (jtwt > max_jtwt) max_jtwt = jtwt;

408:     pre_epsilon = epsilon;
409:     epsilon     = min_jw - jtwt;

411:     if (rank == 0) {
412:       if (innerSolverTol > epsilon) innerSolverTol = epsilon;
413:       else if (innerSolverTol < 1e-7) innerSolverTol = 1e-7;

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

418:       df.tol = innerSolverTol * 0.5;
419:     }

421:     tao->niter++;
422:     PetscCall(TaoLogConvergenceHistory(tao, min_jw, epsilon, 0.0, tao->ksp_its));
423:     PetscCall(TaoMonitor(tao, tao->niter, min_jw, epsilon, 0.0, tao->step));
424:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
425:   }

427:   /* free all the memory */
428:   if (rank == 0) {
429:     PetscCall(destroy_grad_list(&grad_list));
430:     PetscCall(destroy_df_solver(&df));
431:   }

433:   PetscCall(VecDestroy(&bmrm->local_w));
434:   PetscCall(VecScatterDestroy(&bmrm->scatter));
435:   PetscFunctionReturn(PETSC_SUCCESS);
436: }

438: /* ---------------------------------------------------------- */

440: static PetscErrorCode TaoSetup_BMRM(Tao tao)
441: {
442:   PetscFunctionBegin;
443:   /* Allocate some arrays */
444:   if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
445:   PetscFunctionReturn(PETSC_SUCCESS);
446: }

448: /*------------------------------------------------------------*/
449: static PetscErrorCode TaoDestroy_BMRM(Tao tao)
450: {
451:   PetscFunctionBegin;
452:   PetscCall(PetscFree(tao->data));
453:   PetscFunctionReturn(PETSC_SUCCESS);
454: }

456: static PetscErrorCode TaoSetFromOptions_BMRM(Tao tao, PetscOptionItems *PetscOptionsObject)
457: {
458:   TAO_BMRM *bmrm = (TAO_BMRM *)tao->data;

460:   PetscFunctionBegin;
461:   PetscOptionsHeadBegin(PetscOptionsObject, "BMRM for regularized risk minimization");
462:   PetscCall(PetscOptionsReal("-tao_bmrm_lambda", "regulariser weight", "", 100, &bmrm->lambda, NULL));
463:   PetscOptionsHeadEnd();
464:   PetscFunctionReturn(PETSC_SUCCESS);
465: }

467: /*------------------------------------------------------------*/
468: static PetscErrorCode TaoView_BMRM(Tao tao, PetscViewer viewer)
469: {
470:   PetscBool isascii;

472:   PetscFunctionBegin;
473:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
474:   if (isascii) {
475:     PetscCall(PetscViewerASCIIPushTab(viewer));
476:     PetscCall(PetscViewerASCIIPopTab(viewer));
477:   }
478:   PetscFunctionReturn(PETSC_SUCCESS);
479: }

481: /*------------------------------------------------------------*/
482: /*MC
483:   TAOBMRM - bundle method for regularized risk minimization

485:   Options Database Keys:
486: . - tao_bmrm_lambda - regulariser weight

488:   Level: beginner
489: M*/

491: PETSC_EXTERN PetscErrorCode TaoCreate_BMRM(Tao tao)
492: {
493:   TAO_BMRM *bmrm;

495:   PetscFunctionBegin;
496:   tao->ops->setup          = TaoSetup_BMRM;
497:   tao->ops->solve          = TaoSolve_BMRM;
498:   tao->ops->view           = TaoView_BMRM;
499:   tao->ops->setfromoptions = TaoSetFromOptions_BMRM;
500:   tao->ops->destroy        = TaoDestroy_BMRM;

502:   PetscCall(PetscNew(&bmrm));
503:   bmrm->lambda = 1.0;
504:   tao->data    = (void *)bmrm;

506:   /* Override default settings (unless already changed) */
507:   if (!tao->max_it_changed) tao->max_it = 2000;
508:   if (!tao->max_funcs_changed) tao->max_funcs = 4000;
509:   if (!tao->gatol_changed) tao->gatol = 1.0e-12;
510:   if (!tao->grtol_changed) tao->grtol = 1.0e-12;

512:   PetscFunctionReturn(PETSC_SUCCESS);
513: }

515: static PetscErrorCode init_df_solver(TAO_DF *df)
516: {
517:   PetscInt i, n = INCRE_DIM;

519:   PetscFunctionBegin;
520:   /* default values */
521:   df->maxProjIter = 200;
522:   df->maxPGMIter  = 300000;
523:   df->b           = 1.0;

525:   /* memory space required by Dai-Fletcher */
526:   df->cur_num_cp = n;
527:   PetscCall(PetscMalloc1(n, &df->f));
528:   PetscCall(PetscMalloc1(n, &df->a));
529:   PetscCall(PetscMalloc1(n, &df->l));
530:   PetscCall(PetscMalloc1(n, &df->u));
531:   PetscCall(PetscMalloc1(n, &df->x));
532:   PetscCall(PetscMalloc1(n, &df->Q));

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

536:   PetscCall(PetscMalloc1(n, &df->g));
537:   PetscCall(PetscMalloc1(n, &df->y));
538:   PetscCall(PetscMalloc1(n, &df->tempv));
539:   PetscCall(PetscMalloc1(n, &df->d));
540:   PetscCall(PetscMalloc1(n, &df->Qd));
541:   PetscCall(PetscMalloc1(n, &df->t));
542:   PetscCall(PetscMalloc1(n, &df->xplus));
543:   PetscCall(PetscMalloc1(n, &df->tplus));
544:   PetscCall(PetscMalloc1(n, &df->sk));
545:   PetscCall(PetscMalloc1(n, &df->yk));

547:   PetscCall(PetscMalloc1(n, &df->ipt));
548:   PetscCall(PetscMalloc1(n, &df->ipt2));
549:   PetscCall(PetscMalloc1(n, &df->uv));
550:   PetscFunctionReturn(PETSC_SUCCESS);
551: }

553: static PetscErrorCode ensure_df_space(PetscInt dim, TAO_DF *df)
554: {
555:   PetscReal *tmp, **tmp_Q;
556:   PetscInt   i, n, old_n;

558:   PetscFunctionBegin;
559:   df->dim = dim;
560:   if (dim <= df->cur_num_cp) PetscFunctionReturn(PETSC_SUCCESS);

562:   old_n = df->cur_num_cp;
563:   df->cur_num_cp += INCRE_DIM;
564:   n = df->cur_num_cp;

566:   /* memory space required by dai-fletcher */
567:   PetscCall(PetscMalloc1(n, &tmp));
568:   PetscCall(PetscArraycpy(tmp, df->f, old_n));
569:   PetscCall(PetscFree(df->f));
570:   df->f = tmp;

572:   PetscCall(PetscMalloc1(n, &tmp));
573:   PetscCall(PetscArraycpy(tmp, df->a, old_n));
574:   PetscCall(PetscFree(df->a));
575:   df->a = tmp;

577:   PetscCall(PetscMalloc1(n, &tmp));
578:   PetscCall(PetscArraycpy(tmp, df->l, old_n));
579:   PetscCall(PetscFree(df->l));
580:   df->l = tmp;

582:   PetscCall(PetscMalloc1(n, &tmp));
583:   PetscCall(PetscArraycpy(tmp, df->u, old_n));
584:   PetscCall(PetscFree(df->u));
585:   df->u = tmp;

587:   PetscCall(PetscMalloc1(n, &tmp));
588:   PetscCall(PetscArraycpy(tmp, df->x, old_n));
589:   PetscCall(PetscFree(df->x));
590:   df->x = tmp;

592:   PetscCall(PetscMalloc1(n, &tmp_Q));
593:   for (i = 0; i < n; i++) {
594:     PetscCall(PetscMalloc1(n, &tmp_Q[i]));
595:     if (i < old_n) {
596:       PetscCall(PetscArraycpy(tmp_Q[i], df->Q[i], old_n));
597:       PetscCall(PetscFree(df->Q[i]));
598:     }
599:   }

601:   PetscCall(PetscFree(df->Q));
602:   df->Q = tmp_Q;

604:   PetscCall(PetscFree(df->g));
605:   PetscCall(PetscMalloc1(n, &df->g));

607:   PetscCall(PetscFree(df->y));
608:   PetscCall(PetscMalloc1(n, &df->y));

610:   PetscCall(PetscFree(df->tempv));
611:   PetscCall(PetscMalloc1(n, &df->tempv));

613:   PetscCall(PetscFree(df->d));
614:   PetscCall(PetscMalloc1(n, &df->d));

616:   PetscCall(PetscFree(df->Qd));
617:   PetscCall(PetscMalloc1(n, &df->Qd));

619:   PetscCall(PetscFree(df->t));
620:   PetscCall(PetscMalloc1(n, &df->t));

622:   PetscCall(PetscFree(df->xplus));
623:   PetscCall(PetscMalloc1(n, &df->xplus));

625:   PetscCall(PetscFree(df->tplus));
626:   PetscCall(PetscMalloc1(n, &df->tplus));

628:   PetscCall(PetscFree(df->sk));
629:   PetscCall(PetscMalloc1(n, &df->sk));

631:   PetscCall(PetscFree(df->yk));
632:   PetscCall(PetscMalloc1(n, &df->yk));

634:   PetscCall(PetscFree(df->ipt));
635:   PetscCall(PetscMalloc1(n, &df->ipt));

637:   PetscCall(PetscFree(df->ipt2));
638:   PetscCall(PetscMalloc1(n, &df->ipt2));

640:   PetscCall(PetscFree(df->uv));
641:   PetscCall(PetscMalloc1(n, &df->uv));
642:   PetscFunctionReturn(PETSC_SUCCESS);
643: }

645: static PetscErrorCode destroy_df_solver(TAO_DF *df)
646: {
647:   PetscInt i;

649:   PetscFunctionBegin;
650:   PetscCall(PetscFree(df->f));
651:   PetscCall(PetscFree(df->a));
652:   PetscCall(PetscFree(df->l));
653:   PetscCall(PetscFree(df->u));
654:   PetscCall(PetscFree(df->x));

656:   for (i = 0; i < df->cur_num_cp; i++) PetscCall(PetscFree(df->Q[i]));
657:   PetscCall(PetscFree(df->Q));
658:   PetscCall(PetscFree(df->ipt));
659:   PetscCall(PetscFree(df->ipt2));
660:   PetscCall(PetscFree(df->uv));
661:   PetscCall(PetscFree(df->g));
662:   PetscCall(PetscFree(df->y));
663:   PetscCall(PetscFree(df->tempv));
664:   PetscCall(PetscFree(df->d));
665:   PetscCall(PetscFree(df->Qd));
666:   PetscCall(PetscFree(df->t));
667:   PetscCall(PetscFree(df->xplus));
668:   PetscCall(PetscFree(df->tplus));
669:   PetscCall(PetscFree(df->sk));
670:   PetscCall(PetscFree(df->yk));
671:   PetscFunctionReturn(PETSC_SUCCESS);
672: }

674: /* Piecewise linear monotone target function for the Dai-Fletcher projector */
675: static PetscReal phi(PetscReal *x, PetscInt n, PetscReal lambda, PetscReal *a, PetscReal b, PetscReal *c, PetscReal *l, PetscReal *u)
676: {
677:   PetscReal r = 0.0;
678:   PetscInt  i;

680:   for (i = 0; i < n; i++) {
681:     x[i] = -c[i] + lambda * a[i];
682:     if (x[i] > u[i]) x[i] = u[i];
683:     else if (x[i] < l[i]) x[i] = l[i];
684:     r += a[i] * x[i];
685:   }
686:   return r - b;
687: }

689: /** Modified Dai-Fletcher QP projector solves the problem:
690:  *
691:  *      minimise  0.5*x'*x - c'*x
692:  *      subj to   a'*x = b
693:  *                l \leq x \leq u
694:  *
695:  *  \param c The point to be projected onto feasible set
696:  */
697: static PetscInt project(PetscInt n, PetscReal *a, PetscReal b, PetscReal *c, PetscReal *l, PetscReal *u, PetscReal *x, PetscReal *lam_ext, TAO_DF *df)
698: {
699:   PetscReal lambda, lambdal, lambdau, dlambda, lambda_new;
700:   PetscReal r, rl, ru, s;
701:   PetscInt  innerIter;
702:   PetscBool nonNegativeSlack = PETSC_FALSE;

704:   *lam_ext  = 0;
705:   lambda    = 0;
706:   dlambda   = 0.5;
707:   innerIter = 1;

709:   /*  \phi(x;lambda) := 0.5*x'*x + c'*x - lambda*(a'*x-b)
710:    *
711:    *  Optimality conditions for \phi:
712:    *
713:    *  1. lambda   <= 0
714:    *  2. r        <= 0
715:    *  3. r*lambda == 0
716:    */

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

721:   if (nonNegativeSlack) {
722:     /* inequality constraint, i.e., with \xi >= 0 constraint */
723:     if (r < TOL_R) return PETSC_SUCCESS;
724:   } else {
725:     /* equality constraint ,i.e., without \xi >= 0 constraint */
726:     if (PetscAbsReal(r) < TOL_R) return PETSC_SUCCESS;
727:   }

729:   if (r < 0.0) {
730:     lambdal = lambda;
731:     rl      = r;
732:     lambda  = lambda + dlambda;
733:     r       = phi(x, n, lambda, a, b, c, l, u);
734:     while (r < 0.0 && dlambda < BMRM_INFTY) {
735:       lambdal = lambda;
736:       s       = rl / r - 1.0;
737:       if (s < 0.1) s = 0.1;
738:       dlambda = dlambda + dlambda / s;
739:       lambda  = lambda + dlambda;
740:       rl      = r;
741:       r       = phi(x, n, lambda, a, b, c, l, u);
742:     }
743:     lambdau = lambda;
744:     ru      = r;
745:   } else {
746:     lambdau = lambda;
747:     ru      = r;
748:     lambda  = lambda - dlambda;
749:     r       = phi(x, n, lambda, a, b, c, l, u);
750:     while (r > 0.0 && dlambda > -BMRM_INFTY) {
751:       lambdau = lambda;
752:       s       = ru / r - 1.0;
753:       if (s < 0.1) s = 0.1;
754:       dlambda = dlambda + dlambda / s;
755:       lambda  = lambda - dlambda;
756:       ru      = r;
757:       r       = phi(x, n, lambda, a, b, c, l, u);
758:     }
759:     lambdal = lambda;
760:     rl      = r;
761:   }

763:   PetscCheck(PetscAbsReal(dlambda) <= BMRM_INFTY, PETSC_COMM_SELF, PETSC_ERR_PLIB, "L2N2_DaiFletcherPGM detected Infeasible QP problem!");

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

767:   /* Secant Phase */
768:   s       = 1.0 - rl / ru;
769:   dlambda = dlambda / s;
770:   lambda  = lambdau - dlambda;
771:   r       = phi(x, n, lambda, a, b, c, l, u);

773:   while (PetscAbsReal(r) > TOL_R && dlambda > TOL_LAM * (1.0 + PetscAbsReal(lambda)) && innerIter < df->maxProjIter) {
774:     innerIter++;
775:     if (r > 0.0) {
776:       if (s <= 2.0) {
777:         lambdau = lambda;
778:         ru      = r;
779:         s       = 1.0 - rl / ru;
780:         dlambda = (lambdau - lambdal) / s;
781:         lambda  = lambdau - dlambda;
782:       } else {
783:         s = ru / r - 1.0;
784:         if (s < 0.1) s = 0.1;
785:         dlambda    = (lambdau - lambda) / s;
786:         lambda_new = 0.75 * lambdal + 0.25 * lambda;
787:         if (lambda_new < (lambda - dlambda)) lambda_new = lambda - dlambda;
788:         lambdau = lambda;
789:         ru      = r;
790:         lambda  = lambda_new;
791:         s       = (lambdau - lambdal) / (lambdau - lambda);
792:       }
793:     } else {
794:       if (s >= 2.0) {
795:         lambdal = lambda;
796:         rl      = r;
797:         s       = 1.0 - rl / ru;
798:         dlambda = (lambdau - lambdal) / s;
799:         lambda  = lambdau - dlambda;
800:       } else {
801:         s = rl / r - 1.0;
802:         if (s < 0.1) s = 0.1;
803:         dlambda    = (lambda - lambdal) / s;
804:         lambda_new = 0.75 * lambdau + 0.25 * lambda;
805:         if (lambda_new > (lambda + dlambda)) lambda_new = lambda + dlambda;
806:         lambdal = lambda;
807:         rl      = r;
808:         lambda  = lambda_new;
809:         s       = (lambdau - lambdal) / (lambdau - lambda);
810:       }
811:     }
812:     r = phi(x, n, lambda, a, b, c, l, u);
813:   }

815:   *lam_ext = lambda;
816:   if (innerIter >= df->maxProjIter) PetscCall(PetscInfo(NULL, "WARNING: DaiFletcher max iterations\n"));
817:   return innerIter;
818: }