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 PetscCheck(dim >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected dim %" PetscInt_FMT " >= 0", dim); 22 /* variables for the adaptive nonmonotone linesearch */ 23 PetscInt L, llast; 24 PetscReal fr, fbest, fv, fc, fv0; 25 26 c = BMRM_INFTY; 27 28 DELTAsv = EPS_SV; 29 if (tol <= 1.0e-5 || dim <= 20) ProdDELTAsv = 0.0F; 30 else ProdDELTAsv = EPS_SV; 31 32 for (i = 0; i < dim; i++) tempv[i] = -x[i]; 33 34 lam_ext = 0.0; 35 36 /* Project the initial solution */ 37 project(dim, a, b, tempv, l, u, x, &lam_ext, df); 38 39 /* Compute gradient 40 g = Q*x + f; */ 41 42 it = 0; 43 for (i = 0; i < dim; i++) { 44 if (PetscAbsReal(x[i]) > ProdDELTAsv) ipt[it++] = i; 45 } 46 47 PetscCall(PetscArrayzero(t, dim)); 48 for (i = 0; i < it; i++) { 49 tempQ = Q[ipt[i]]; 50 for (j = 0; j < dim; j++) t[j] += (tempQ[j] * x[ipt[i]]); 51 } 52 for (i = 0; i < dim; i++) g[i] = t[i] + f[i]; 53 54 /* y = -(x_{k} - g_{k}) */ 55 for (i = 0; i < dim; i++) y[i] = g[i] - x[i]; 56 57 /* Project x_{k} - g_{k} */ 58 project(dim, a, b, y, l, u, tempv, &lam_ext, df); 59 60 /* y = P(x_{k} - g_{k}) - x_{k} */ 61 max = ALPHA_MIN; 62 for (i = 0; i < dim; i++) { 63 y[i] = tempv[i] - x[i]; 64 if (PetscAbsReal(y[i]) > max) max = PetscAbsReal(y[i]); 65 } 66 67 if (max < tol * 1e-3) return PETSC_SUCCESS; 68 69 alpha = 1.0 / max; 70 71 /* fv0 = f(x_{0}). Recall t = Q x_{k} */ 72 fv0 = 0.0; 73 for (i = 0; i < dim; i++) fv0 += x[i] * (0.5 * t[i] + f[i]); 74 75 /* adaptive nonmonotone linesearch */ 76 L = 2; 77 fr = ALPHA_MAX; 78 fbest = fv0; 79 fc = fv0; 80 llast = 0; 81 akold = bkold = 0.0; 82 83 /* Iterator begins */ 84 for (innerIter = 1; innerIter <= df->maxPGMIter; innerIter++) { 85 /* tempv = -(x_{k} - alpha*g_{k}) */ 86 for (i = 0; i < dim; i++) tempv[i] = alpha * g[i] - x[i]; 87 88 /* Project x_{k} - alpha*g_{k} */ 89 project(dim, a, b, tempv, l, u, y, &lam_ext, df); 90 91 /* gd = \inner{d_{k}}{g_{k}} 92 d = P(x_{k} - alpha*g_{k}) - x_{k} 93 */ 94 gd = 0.0; 95 for (i = 0; i < dim; i++) { 96 d[i] = y[i] - x[i]; 97 gd += d[i] * g[i]; 98 } 99 100 /* Gradient computation */ 101 102 /* compute Qd = Q*d or Qd = Q*y - t depending on their sparsity */ 103 104 it = it2 = 0; 105 for (i = 0; i < dim; i++) { 106 if (PetscAbsReal(d[i]) > (ProdDELTAsv * 1.0e-2)) ipt[it++] = i; 107 } 108 for (i = 0; i < dim; i++) { 109 if (PetscAbsReal(y[i]) > ProdDELTAsv) ipt2[it2++] = i; 110 } 111 112 PetscCall(PetscArrayzero(Qd, dim)); 113 /* compute Qd = Q*d */ 114 if (it < it2) { 115 for (i = 0; i < it; i++) { 116 tempQ = Q[ipt[i]]; 117 for (j = 0; j < dim; j++) Qd[j] += (tempQ[j] * d[ipt[i]]); 118 } 119 } else { /* compute Qd = Q*y-t */ 120 for (i = 0; i < it2; i++) { 121 tempQ = Q[ipt2[i]]; 122 for (j = 0; j < dim; j++) Qd[j] += (tempQ[j] * y[ipt2[i]]); 123 } 124 for (j = 0; j < dim; j++) Qd[j] -= t[j]; 125 } 126 127 /* ak = inner{d_{k}}{d_{k}} */ 128 ak = 0.0; 129 for (i = 0; i < dim; i++) ak += d[i] * d[i]; 130 131 bk = 0.0; 132 for (i = 0; i < dim; i++) bk += d[i] * Qd[i]; 133 134 if (bk > EPS * ak && gd < 0.0) lamnew = -gd / bk; 135 else lamnew = 1.0; 136 137 /* fv is computing f(x_{k} + d_{k}) */ 138 fv = 0.0; 139 for (i = 0; i < dim; i++) { 140 xplus[i] = x[i] + d[i]; 141 tplus[i] = t[i] + Qd[i]; 142 fv += xplus[i] * (0.5 * tplus[i] + f[i]); 143 } 144 145 /* fr is fref */ 146 if ((innerIter == 1 && fv >= fv0) || (innerIter > 1 && fv >= fr)) { 147 fv = 0.0; 148 for (i = 0; i < dim; i++) { 149 xplus[i] = x[i] + lamnew * d[i]; 150 tplus[i] = t[i] + lamnew * Qd[i]; 151 fv += xplus[i] * (0.5 * tplus[i] + f[i]); 152 } 153 } 154 155 for (i = 0; i < dim; i++) { 156 sk[i] = xplus[i] - x[i]; 157 yk[i] = tplus[i] - t[i]; 158 x[i] = xplus[i]; 159 t[i] = tplus[i]; 160 g[i] = t[i] + f[i]; 161 } 162 163 /* update the line search control parameters */ 164 if (fv < fbest) { 165 fbest = fv; 166 fc = fv; 167 llast = 0; 168 } else { 169 fc = (fc > fv ? fc : fv); 170 llast++; 171 if (llast == L) { 172 fr = fc; 173 fc = fv; 174 llast = 0; 175 } 176 } 177 178 ak = bk = 0.0; 179 for (i = 0; i < dim; i++) { 180 ak += sk[i] * sk[i]; 181 bk += sk[i] * yk[i]; 182 } 183 184 if (bk <= EPS * ak) alpha = ALPHA_MAX; 185 else { 186 if (bkold < EPS * akold) alpha = ak / bk; 187 else alpha = (akold + ak) / (bkold + bk); 188 189 if (alpha > ALPHA_MAX) alpha = ALPHA_MAX; 190 else if (alpha < ALPHA_MIN) alpha = ALPHA_MIN; 191 } 192 193 akold = ak; 194 bkold = bk; 195 196 /* stopping criterion based on KKT conditions */ 197 /* at optimal, gradient of lagrangian w.r.t. x is zero */ 198 199 bk = 0.0; 200 for (i = 0; i < dim; i++) bk += x[i] * x[i]; 201 202 if (PetscSqrtReal(ak) < tol * 10 * PetscSqrtReal(bk)) { 203 it = 0; 204 luv = 0; 205 kktlam = 0.0; 206 for (i = 0; i < dim; i++) { 207 /* x[i] is active hence lagrange multipliers for box constraints 208 are zero. The lagrange multiplier for ineq. const. is then 209 defined as below 210 */ 211 if ((x[i] > DELTAsv) && (x[i] < c - DELTAsv)) { 212 ipt[it++] = i; 213 kktlam = kktlam - a[i] * g[i]; 214 } else uv[luv++] = i; 215 } 216 217 if (it == 0 && PetscSqrtReal(ak) < tol * 0.5 * PetscSqrtReal(bk)) return PETSC_SUCCESS; 218 else { 219 kktlam = kktlam / it; 220 info = 1; 221 for (i = 0; i < it; i++) { 222 if (PetscAbsReal(a[ipt[i]] * g[ipt[i]] + kktlam) > tol) { 223 info = 0; 224 break; 225 } 226 } 227 if (info == 1) { 228 for (i = 0; i < luv; i++) { 229 if (x[uv[i]] <= DELTAsv) { 230 /* x[i] == lower bound, hence, lagrange multiplier (say, beta) for lower bound may 231 not be zero. So, the gradient without beta is > 0 232 */ 233 if (g[uv[i]] + kktlam * a[uv[i]] < -tol) { 234 info = 0; 235 break; 236 } 237 } else { 238 /* x[i] == upper bound, hence, lagrange multiplier (say, eta) for upper bound may 239 not be zero. So, the gradient without eta is < 0 240 */ 241 if (g[uv[i]] + kktlam * a[uv[i]] > tol) { 242 info = 0; 243 break; 244 } 245 } 246 } 247 } 248 249 if (info == 1) return PETSC_SUCCESS; 250 } 251 } 252 } 253 return PETSC_SUCCESS; 254 } 255 256 /* The main solver function 257 258 f = Remp(W) This is what the user provides us from the application layer 259 So the ComputeGradient function for instance should get us back the subgradient of Remp(W) 260 261 Regularizer assumed to be L2 norm = lambda*0.5*W'W () 262 */ 263 264 static PetscErrorCode make_grad_node(Vec X, Vec_Chain **p) 265 { 266 PetscFunctionBegin; 267 PetscCall(PetscNew(p)); 268 PetscCall(VecDuplicate(X, &(*p)->V)); 269 PetscCall(VecCopy(X, (*p)->V)); 270 (*p)->next = NULL; 271 PetscFunctionReturn(PETSC_SUCCESS); 272 } 273 274 static PetscErrorCode destroy_grad_list(Vec_Chain *head) 275 { 276 Vec_Chain *p = head->next, *q; 277 278 PetscFunctionBegin; 279 while (p) { 280 q = p->next; 281 PetscCall(VecDestroy(&p->V)); 282 PetscCall(PetscFree(p)); 283 p = q; 284 } 285 head->next = NULL; 286 PetscFunctionReturn(PETSC_SUCCESS); 287 } 288 289 static PetscErrorCode TaoSolve_BMRM(Tao tao) 290 { 291 TAO_DF df; 292 TAO_BMRM *bmrm = (TAO_BMRM *)tao->data; 293 294 /* Values and pointers to parts of the optimization problem */ 295 PetscReal f = 0.0; 296 Vec W = tao->solution; 297 Vec G = tao->gradient; 298 PetscReal lambda; 299 PetscReal bt; 300 Vec_Chain grad_list, *tail_glist, *pgrad; 301 PetscInt i; 302 PetscMPIInt rank; 303 304 /* Used in converged criteria check */ 305 PetscReal reg; 306 PetscReal jtwt = 0.0, max_jtwt, pre_epsilon, epsilon, jw, min_jw; 307 PetscReal innerSolverTol; 308 MPI_Comm comm; 309 310 PetscFunctionBegin; 311 PetscCall(PetscObjectGetComm((PetscObject)tao, &comm)); 312 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 313 lambda = bmrm->lambda; 314 315 /* Check Stopping Condition */ 316 tao->step = 1.0; 317 max_jtwt = -BMRM_INFTY; 318 min_jw = BMRM_INFTY; 319 innerSolverTol = 1.0; 320 epsilon = 0.0; 321 322 if (rank == 0) { 323 PetscCall(init_df_solver(&df)); 324 grad_list.next = NULL; 325 tail_glist = &grad_list; 326 } 327 328 df.tol = 1e-6; 329 tao->reason = TAO_CONTINUE_ITERATING; 330 331 /*-----------------Algorithm Begins------------------------*/ 332 /* make the scatter */ 333 PetscCall(VecScatterCreateToZero(W, &bmrm->scatter, &bmrm->local_w)); 334 PetscCall(VecAssemblyBegin(bmrm->local_w)); 335 PetscCall(VecAssemblyEnd(bmrm->local_w)); 336 337 /* NOTE: In application pass the sub-gradient of Remp(W) */ 338 PetscCall(TaoComputeObjectiveAndGradient(tao, W, &f, G)); 339 PetscCall(TaoLogConvergenceHistory(tao, f, 1.0, 0.0, tao->ksp_its)); 340 PetscCall(TaoMonitor(tao, tao->niter, f, 1.0, 0.0, tao->step)); 341 PetscUseTypeMethod(tao, convergencetest, tao->cnvP); 342 343 while (tao->reason == TAO_CONTINUE_ITERATING) { 344 /* Call general purpose update function */ 345 PetscTryTypeMethod(tao, update, tao->niter, tao->user_update); 346 347 /* compute bt = Remp(Wt-1) - <Wt-1, At> */ 348 PetscCall(VecDot(W, G, &bt)); 349 bt = f - bt; 350 351 /* First gather the gradient to the rank-0 node */ 352 PetscCall(VecScatterBegin(bmrm->scatter, G, bmrm->local_w, INSERT_VALUES, SCATTER_FORWARD)); 353 PetscCall(VecScatterEnd(bmrm->scatter, G, bmrm->local_w, INSERT_VALUES, SCATTER_FORWARD)); 354 355 /* Bring up the inner solver */ 356 if (rank == 0) { 357 PetscCall(ensure_df_space(tao->niter + 1, &df)); 358 PetscCall(make_grad_node(bmrm->local_w, &pgrad)); 359 tail_glist->next = pgrad; 360 tail_glist = pgrad; 361 362 df.a[tao->niter] = 1.0; 363 df.f[tao->niter] = -bt; 364 df.u[tao->niter] = 1.0; 365 df.l[tao->niter] = 0.0; 366 367 /* set up the Q */ 368 pgrad = grad_list.next; 369 for (i = 0; i <= tao->niter; i++) { 370 PetscCheck(pgrad, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Assert that there are at least tao->niter+1 pgrad available"); 371 PetscCall(VecDot(pgrad->V, bmrm->local_w, ®)); 372 df.Q[i][tao->niter] = df.Q[tao->niter][i] = reg / lambda; 373 pgrad = pgrad->next; 374 } 375 376 if (tao->niter > 0) { 377 df.x[tao->niter] = 0.0; 378 PetscCall(solve(&df)); 379 } else df.x[0] = 1.0; 380 381 /* now computing Jt*(alpha_t) which should be = Jt(wt) to check convergence */ 382 jtwt = 0.0; 383 PetscCall(VecSet(bmrm->local_w, 0.0)); 384 pgrad = grad_list.next; 385 for (i = 0; i <= tao->niter; i++) { 386 jtwt -= df.x[i] * df.f[i]; 387 PetscCall(VecAXPY(bmrm->local_w, -df.x[i] / lambda, pgrad->V)); 388 pgrad = pgrad->next; 389 } 390 391 PetscCall(VecNorm(bmrm->local_w, NORM_2, ®)); 392 reg = 0.5 * lambda * reg * reg; 393 jtwt -= reg; 394 } /* end if rank == 0 */ 395 396 /* scatter the new W to all nodes */ 397 PetscCall(VecScatterBegin(bmrm->scatter, bmrm->local_w, W, INSERT_VALUES, SCATTER_REVERSE)); 398 PetscCall(VecScatterEnd(bmrm->scatter, bmrm->local_w, W, INSERT_VALUES, SCATTER_REVERSE)); 399 400 PetscCall(TaoComputeObjectiveAndGradient(tao, W, &f, G)); 401 402 PetscCallMPI(MPI_Bcast(&jtwt, 1, MPIU_REAL, 0, comm)); 403 PetscCallMPI(MPI_Bcast(®, 1, MPIU_REAL, 0, comm)); 404 405 jw = reg + f; /* J(w) = regularizer + Remp(w) */ 406 if (jw < min_jw) min_jw = jw; 407 if (jtwt > max_jtwt) max_jtwt = jtwt; 408 409 pre_epsilon = epsilon; 410 epsilon = min_jw - jtwt; 411 412 if (rank == 0) { 413 if (innerSolverTol > epsilon) innerSolverTol = epsilon; 414 else if (innerSolverTol < 1e-7) innerSolverTol = 1e-7; 415 416 /* if the annealing doesn't work well, lower the inner solver tolerance */ 417 if (pre_epsilon < epsilon) innerSolverTol *= 0.2; 418 419 df.tol = innerSolverTol * 0.5; 420 } 421 422 tao->niter++; 423 PetscCall(TaoLogConvergenceHistory(tao, min_jw, epsilon, 0.0, tao->ksp_its)); 424 PetscCall(TaoMonitor(tao, tao->niter, min_jw, epsilon, 0.0, tao->step)); 425 PetscUseTypeMethod(tao, convergencetest, tao->cnvP); 426 } 427 428 /* free all the memory */ 429 if (rank == 0) { 430 PetscCall(destroy_grad_list(&grad_list)); 431 PetscCall(destroy_df_solver(&df)); 432 } 433 434 PetscCall(VecDestroy(&bmrm->local_w)); 435 PetscCall(VecScatterDestroy(&bmrm->scatter)); 436 PetscFunctionReturn(PETSC_SUCCESS); 437 } 438 439 /* ---------------------------------------------------------- */ 440 441 static PetscErrorCode TaoSetup_BMRM(Tao tao) 442 { 443 PetscFunctionBegin; 444 /* Allocate some arrays */ 445 if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient)); 446 PetscFunctionReturn(PETSC_SUCCESS); 447 } 448 449 /*------------------------------------------------------------*/ 450 static PetscErrorCode TaoDestroy_BMRM(Tao tao) 451 { 452 PetscFunctionBegin; 453 PetscCall(PetscFree(tao->data)); 454 PetscFunctionReturn(PETSC_SUCCESS); 455 } 456 457 static PetscErrorCode TaoSetFromOptions_BMRM(Tao tao, PetscOptionItems *PetscOptionsObject) 458 { 459 TAO_BMRM *bmrm = (TAO_BMRM *)tao->data; 460 461 PetscFunctionBegin; 462 PetscOptionsHeadBegin(PetscOptionsObject, "BMRM for regularized risk minimization"); 463 PetscCall(PetscOptionsReal("-tao_bmrm_lambda", "regulariser weight", "", 100, &bmrm->lambda, NULL)); 464 PetscOptionsHeadEnd(); 465 PetscFunctionReturn(PETSC_SUCCESS); 466 } 467 468 /*------------------------------------------------------------*/ 469 static PetscErrorCode TaoView_BMRM(Tao tao, PetscViewer viewer) 470 { 471 PetscBool isascii; 472 473 PetscFunctionBegin; 474 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 475 if (isascii) { 476 PetscCall(PetscViewerASCIIPushTab(viewer)); 477 PetscCall(PetscViewerASCIIPopTab(viewer)); 478 } 479 PetscFunctionReturn(PETSC_SUCCESS); 480 } 481 482 /*------------------------------------------------------------*/ 483 /*MC 484 TAOBMRM - bundle method for regularized risk minimization 485 486 Options Database Keys: 487 . - tao_bmrm_lambda - regulariser weight 488 489 Level: beginner 490 M*/ 491 492 PETSC_EXTERN PetscErrorCode TaoCreate_BMRM(Tao tao) 493 { 494 TAO_BMRM *bmrm; 495 496 PetscFunctionBegin; 497 tao->ops->setup = TaoSetup_BMRM; 498 tao->ops->solve = TaoSolve_BMRM; 499 tao->ops->view = TaoView_BMRM; 500 tao->ops->setfromoptions = TaoSetFromOptions_BMRM; 501 tao->ops->destroy = TaoDestroy_BMRM; 502 503 PetscCall(PetscNew(&bmrm)); 504 bmrm->lambda = 1.0; 505 tao->data = (void *)bmrm; 506 507 /* Override default settings (unless already changed) */ 508 if (!tao->max_it_changed) tao->max_it = 2000; 509 if (!tao->max_funcs_changed) tao->max_funcs = 4000; 510 if (!tao->gatol_changed) tao->gatol = 1.0e-12; 511 if (!tao->grtol_changed) tao->grtol = 1.0e-12; 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