1 #include <../src/tao/unconstrained/impls/bmrm/bmrm.h> 2 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 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; 20 21 /* variables for the adaptive nonmonotone linesearch */ 22 PetscInt L, llast; 23 PetscReal fr, fbest, fv, fc, fv0; 24 25 c = BMRM_INFTY; 26 27 DELTAsv = EPS_SV; 28 if (tol <= 1.0e-5 || dim <= 20) ProdDELTAsv = 0.0F; 29 else ProdDELTAsv = EPS_SV; 30 31 for (i = 0; i < dim; i++) tempv[i] = -x[i]; 32 33 lam_ext = 0.0; 34 35 /* Project the initial solution */ 36 project(dim, a, b, tempv, l, u, x, &lam_ext, df); 37 38 /* Compute gradient 39 g = Q*x + f; */ 40 41 it = 0; 42 for (i = 0; i < dim; i++) { 43 if (PetscAbsReal(x[i]) > ProdDELTAsv) ipt[it++] = i; 44 } 45 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]; 52 53 /* y = -(x_{k} - g_{k}) */ 54 for (i = 0; i < dim; i++) y[i] = g[i] - x[i]; 55 56 /* Project x_{k} - g_{k} */ 57 project(dim, a, b, y, l, u, tempv, &lam_ext, df); 58 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 } 65 66 if (max < tol * 1e-3) return PETSC_SUCCESS; 67 68 alpha = 1.0 / max; 69 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]); 73 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; 81 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]; 86 87 /* Project x_{k} - alpha*g_{k} */ 88 project(dim, a, b, tempv, l, u, y, &lam_ext, df); 89 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 } 98 99 /* Gradient computation */ 100 101 /* compute Qd = Q*d or Qd = Q*y - t depending on their sparsity */ 102 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 } 110 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 } 125 126 /* ak = inner{d_{k}}{d_{k}} */ 127 ak = 0.0; 128 for (i = 0; i < dim; i++) ak += d[i] * d[i]; 129 130 bk = 0.0; 131 for (i = 0; i < dim; i++) bk += d[i] * Qd[i]; 132 133 if (bk > EPS * ak && gd < 0.0) lamnew = -gd / bk; 134 else lamnew = 1.0; 135 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 } 143 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 } 153 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 } 161 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 } 176 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 } 182 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); 187 188 if (alpha > ALPHA_MAX) alpha = ALPHA_MAX; 189 else if (alpha < ALPHA_MIN) alpha = ALPHA_MIN; 190 } 191 192 akold = ak; 193 bkold = bk; 194 195 /* stopping criterion based on KKT conditions */ 196 /* at optimal, gradient of lagrangian w.r.t. x is zero */ 197 198 bk = 0.0; 199 for (i = 0; i < dim; i++) bk += x[i] * x[i]; 200 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 } 215 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 } 247 248 if (info == 1) return PETSC_SUCCESS; 249 } 250 } 251 } 252 return PETSC_SUCCESS; 253 } 254 255 /* The main solver function 256 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) 259 260 Regularizer assumed to be L2 norm = lambda*0.5*W'W () 261 */ 262 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 } 272 273 static PetscErrorCode destroy_grad_list(Vec_Chain *head) 274 { 275 Vec_Chain *p = head->next, *q; 276 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 } 287 288 static PetscErrorCode TaoSolve_BMRM(Tao tao) 289 { 290 TAO_DF df; 291 TAO_BMRM *bmrm = (TAO_BMRM *)tao->data; 292 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; 302 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; 308 309 PetscFunctionBegin; 310 PetscCall(PetscObjectGetComm((PetscObject)tao, &comm)); 311 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 312 lambda = bmrm->lambda; 313 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; 320 321 if (rank == 0) { 322 PetscCall(init_df_solver(&df)); 323 grad_list.next = NULL; 324 tail_glist = &grad_list; 325 } 326 327 df.tol = 1e-6; 328 tao->reason = TAO_CONTINUE_ITERATING; 329 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)); 335 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); 341 342 while (tao->reason == TAO_CONTINUE_ITERATING) { 343 /* Call general purpose update function */ 344 PetscTryTypeMethod(tao, update, tao->niter, tao->user_update); 345 346 /* compute bt = Remp(Wt-1) - <Wt-1, At> */ 347 PetscCall(VecDot(W, G, &bt)); 348 bt = f - bt; 349 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)); 353 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; 360 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; 365 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, ®)); 371 df.Q[i][tao->niter] = df.Q[tao->niter][i] = reg / lambda; 372 pgrad = pgrad->next; 373 } 374 375 if (tao->niter > 0) { 376 df.x[tao->niter] = 0.0; 377 PetscCall(solve(&df)); 378 } else df.x[0] = 1.0; 379 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 } 389 390 PetscCall(VecNorm(bmrm->local_w, NORM_2, ®)); 391 reg = 0.5 * lambda * reg * reg; 392 jtwt -= reg; 393 } /* end if rank == 0 */ 394 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)); 398 399 PetscCall(TaoComputeObjectiveAndGradient(tao, W, &f, G)); 400 401 PetscCallMPI(MPI_Bcast(&jtwt, 1, MPIU_REAL, 0, comm)); 402 PetscCallMPI(MPI_Bcast(®, 1, MPIU_REAL, 0, comm)); 403 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; 407 408 pre_epsilon = epsilon; 409 epsilon = min_jw - jtwt; 410 411 if (rank == 0) { 412 if (innerSolverTol > epsilon) innerSolverTol = epsilon; 413 else if (innerSolverTol < 1e-7) innerSolverTol = 1e-7; 414 415 /* if the annealing doesn't work well, lower the inner solver tolerance */ 416 if (pre_epsilon < epsilon) innerSolverTol *= 0.2; 417 418 df.tol = innerSolverTol * 0.5; 419 } 420 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 } 426 427 /* free all the memory */ 428 if (rank == 0) { 429 PetscCall(destroy_grad_list(&grad_list)); 430 PetscCall(destroy_df_solver(&df)); 431 } 432 433 PetscCall(VecDestroy(&bmrm->local_w)); 434 PetscCall(VecScatterDestroy(&bmrm->scatter)); 435 PetscFunctionReturn(PETSC_SUCCESS); 436 } 437 438 /* ---------------------------------------------------------- */ 439 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 } 447 448 /*------------------------------------------------------------*/ 449 static PetscErrorCode TaoDestroy_BMRM(Tao tao) 450 { 451 PetscFunctionBegin; 452 PetscCall(PetscFree(tao->data)); 453 PetscFunctionReturn(PETSC_SUCCESS); 454 } 455 456 static PetscErrorCode TaoSetFromOptions_BMRM(Tao tao, PetscOptionItems *PetscOptionsObject) 457 { 458 TAO_BMRM *bmrm = (TAO_BMRM *)tao->data; 459 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 } 466 467 /*------------------------------------------------------------*/ 468 static PetscErrorCode TaoView_BMRM(Tao tao, PetscViewer viewer) 469 { 470 PetscBool isascii; 471 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 } 480 481 /*------------------------------------------------------------*/ 482 /*MC 483 TAOBMRM - bundle method for regularized risk minimization 484 485 Options Database Keys: 486 . - tao_bmrm_lambda - regulariser weight 487 488 Level: beginner 489 M*/ 490 491 PETSC_EXTERN PetscErrorCode TaoCreate_BMRM(Tao tao) 492 { 493 TAO_BMRM *bmrm; 494 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; 501 502 PetscCall(PetscNew(&bmrm)); 503 bmrm->lambda = 1.0; 504 tao->data = (void *)bmrm; 505 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; 511 512 PetscFunctionReturn(PETSC_SUCCESS); 513 } 514 515 static PetscErrorCode init_df_solver(TAO_DF *df) 516 { 517 PetscInt i, n = INCRE_DIM; 518 519 PetscFunctionBegin; 520 /* default values */ 521 df->maxProjIter = 200; 522 df->maxPGMIter = 300000; 523 df->b = 1.0; 524 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)); 533 534 for (i = 0; i < n; i++) PetscCall(PetscMalloc1(n, &df->Q[i])); 535 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)); 546 547 PetscCall(PetscMalloc1(n, &df->ipt)); 548 PetscCall(PetscMalloc1(n, &df->ipt2)); 549 PetscCall(PetscMalloc1(n, &df->uv)); 550 PetscFunctionReturn(PETSC_SUCCESS); 551 } 552 553 static PetscErrorCode ensure_df_space(PetscInt dim, TAO_DF *df) 554 { 555 PetscReal *tmp, **tmp_Q; 556 PetscInt i, n, old_n; 557 558 PetscFunctionBegin; 559 df->dim = dim; 560 if (dim <= df->cur_num_cp) PetscFunctionReturn(PETSC_SUCCESS); 561 562 old_n = df->cur_num_cp; 563 df->cur_num_cp += INCRE_DIM; 564 n = df->cur_num_cp; 565 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; 571 572 PetscCall(PetscMalloc1(n, &tmp)); 573 PetscCall(PetscArraycpy(tmp, df->a, old_n)); 574 PetscCall(PetscFree(df->a)); 575 df->a = tmp; 576 577 PetscCall(PetscMalloc1(n, &tmp)); 578 PetscCall(PetscArraycpy(tmp, df->l, old_n)); 579 PetscCall(PetscFree(df->l)); 580 df->l = tmp; 581 582 PetscCall(PetscMalloc1(n, &tmp)); 583 PetscCall(PetscArraycpy(tmp, df->u, old_n)); 584 PetscCall(PetscFree(df->u)); 585 df->u = tmp; 586 587 PetscCall(PetscMalloc1(n, &tmp)); 588 PetscCall(PetscArraycpy(tmp, df->x, old_n)); 589 PetscCall(PetscFree(df->x)); 590 df->x = tmp; 591 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 } 600 601 PetscCall(PetscFree(df->Q)); 602 df->Q = tmp_Q; 603 604 PetscCall(PetscFree(df->g)); 605 PetscCall(PetscMalloc1(n, &df->g)); 606 607 PetscCall(PetscFree(df->y)); 608 PetscCall(PetscMalloc1(n, &df->y)); 609 610 PetscCall(PetscFree(df->tempv)); 611 PetscCall(PetscMalloc1(n, &df->tempv)); 612 613 PetscCall(PetscFree(df->d)); 614 PetscCall(PetscMalloc1(n, &df->d)); 615 616 PetscCall(PetscFree(df->Qd)); 617 PetscCall(PetscMalloc1(n, &df->Qd)); 618 619 PetscCall(PetscFree(df->t)); 620 PetscCall(PetscMalloc1(n, &df->t)); 621 622 PetscCall(PetscFree(df->xplus)); 623 PetscCall(PetscMalloc1(n, &df->xplus)); 624 625 PetscCall(PetscFree(df->tplus)); 626 PetscCall(PetscMalloc1(n, &df->tplus)); 627 628 PetscCall(PetscFree(df->sk)); 629 PetscCall(PetscMalloc1(n, &df->sk)); 630 631 PetscCall(PetscFree(df->yk)); 632 PetscCall(PetscMalloc1(n, &df->yk)); 633 634 PetscCall(PetscFree(df->ipt)); 635 PetscCall(PetscMalloc1(n, &df->ipt)); 636 637 PetscCall(PetscFree(df->ipt2)); 638 PetscCall(PetscMalloc1(n, &df->ipt2)); 639 640 PetscCall(PetscFree(df->uv)); 641 PetscCall(PetscMalloc1(n, &df->uv)); 642 PetscFunctionReturn(PETSC_SUCCESS); 643 } 644 645 static PetscErrorCode destroy_df_solver(TAO_DF *df) 646 { 647 PetscInt i; 648 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)); 655 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 } 673 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; 679 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 } 688 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; 703 704 *lam_ext = 0; 705 lambda = 0; 706 dlambda = 0.5; 707 innerIter = 1; 708 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 */ 717 718 /* Bracketing Phase */ 719 r = phi(x, n, lambda, a, b, c, l, u); 720 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 } 728 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 } 762 763 PetscCheck(PetscAbsReal(dlambda) <= BMRM_INFTY, PETSC_COMM_SELF, PETSC_ERR_PLIB, "L2N2_DaiFletcherPGM detected Infeasible QP problem!"); 764 765 if (ru == 0) return innerIter; 766 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); 772 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 } 814 815 *lam_ext = lambda; 816 if (innerIter >= df->maxProjIter) PetscCall(PetscInfo(NULL, "WARNING: DaiFletcher max iterations\n")); 817 return innerIter; 818 } 819