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 PetscFunctionReturn(PETSC_SUCCESS); 512 } 513 514 static PetscErrorCode init_df_solver(TAO_DF *df) 515 { 516 PetscInt i, n = INCRE_DIM; 517 518 PetscFunctionBegin; 519 /* default values */ 520 df->maxProjIter = 200; 521 df->maxPGMIter = 300000; 522 df->b = 1.0; 523 524 /* memory space required by Dai-Fletcher */ 525 df->cur_num_cp = n; 526 PetscCall(PetscMalloc1(n, &df->f)); 527 PetscCall(PetscMalloc1(n, &df->a)); 528 PetscCall(PetscMalloc1(n, &df->l)); 529 PetscCall(PetscMalloc1(n, &df->u)); 530 PetscCall(PetscMalloc1(n, &df->x)); 531 PetscCall(PetscMalloc1(n, &df->Q)); 532 533 for (i = 0; i < n; i++) PetscCall(PetscMalloc1(n, &df->Q[i])); 534 535 PetscCall(PetscMalloc1(n, &df->g)); 536 PetscCall(PetscMalloc1(n, &df->y)); 537 PetscCall(PetscMalloc1(n, &df->tempv)); 538 PetscCall(PetscMalloc1(n, &df->d)); 539 PetscCall(PetscMalloc1(n, &df->Qd)); 540 PetscCall(PetscMalloc1(n, &df->t)); 541 PetscCall(PetscMalloc1(n, &df->xplus)); 542 PetscCall(PetscMalloc1(n, &df->tplus)); 543 PetscCall(PetscMalloc1(n, &df->sk)); 544 PetscCall(PetscMalloc1(n, &df->yk)); 545 546 PetscCall(PetscMalloc1(n, &df->ipt)); 547 PetscCall(PetscMalloc1(n, &df->ipt2)); 548 PetscCall(PetscMalloc1(n, &df->uv)); 549 PetscFunctionReturn(PETSC_SUCCESS); 550 } 551 552 static PetscErrorCode ensure_df_space(PetscInt dim, TAO_DF *df) 553 { 554 PetscReal *tmp, **tmp_Q; 555 PetscInt i, n, old_n; 556 557 PetscFunctionBegin; 558 df->dim = dim; 559 if (dim <= df->cur_num_cp) PetscFunctionReturn(PETSC_SUCCESS); 560 561 old_n = df->cur_num_cp; 562 df->cur_num_cp += INCRE_DIM; 563 n = df->cur_num_cp; 564 565 /* memory space required by dai-fletcher */ 566 PetscCall(PetscMalloc1(n, &tmp)); 567 PetscCall(PetscArraycpy(tmp, df->f, old_n)); 568 PetscCall(PetscFree(df->f)); 569 df->f = tmp; 570 571 PetscCall(PetscMalloc1(n, &tmp)); 572 PetscCall(PetscArraycpy(tmp, df->a, old_n)); 573 PetscCall(PetscFree(df->a)); 574 df->a = tmp; 575 576 PetscCall(PetscMalloc1(n, &tmp)); 577 PetscCall(PetscArraycpy(tmp, df->l, old_n)); 578 PetscCall(PetscFree(df->l)); 579 df->l = tmp; 580 581 PetscCall(PetscMalloc1(n, &tmp)); 582 PetscCall(PetscArraycpy(tmp, df->u, old_n)); 583 PetscCall(PetscFree(df->u)); 584 df->u = tmp; 585 586 PetscCall(PetscMalloc1(n, &tmp)); 587 PetscCall(PetscArraycpy(tmp, df->x, old_n)); 588 PetscCall(PetscFree(df->x)); 589 df->x = tmp; 590 591 PetscCall(PetscMalloc1(n, &tmp_Q)); 592 for (i = 0; i < n; i++) { 593 PetscCall(PetscMalloc1(n, &tmp_Q[i])); 594 if (i < old_n) { 595 PetscCall(PetscArraycpy(tmp_Q[i], df->Q[i], old_n)); 596 PetscCall(PetscFree(df->Q[i])); 597 } 598 } 599 600 PetscCall(PetscFree(df->Q)); 601 df->Q = tmp_Q; 602 603 PetscCall(PetscFree(df->g)); 604 PetscCall(PetscMalloc1(n, &df->g)); 605 606 PetscCall(PetscFree(df->y)); 607 PetscCall(PetscMalloc1(n, &df->y)); 608 609 PetscCall(PetscFree(df->tempv)); 610 PetscCall(PetscMalloc1(n, &df->tempv)); 611 612 PetscCall(PetscFree(df->d)); 613 PetscCall(PetscMalloc1(n, &df->d)); 614 615 PetscCall(PetscFree(df->Qd)); 616 PetscCall(PetscMalloc1(n, &df->Qd)); 617 618 PetscCall(PetscFree(df->t)); 619 PetscCall(PetscMalloc1(n, &df->t)); 620 621 PetscCall(PetscFree(df->xplus)); 622 PetscCall(PetscMalloc1(n, &df->xplus)); 623 624 PetscCall(PetscFree(df->tplus)); 625 PetscCall(PetscMalloc1(n, &df->tplus)); 626 627 PetscCall(PetscFree(df->sk)); 628 PetscCall(PetscMalloc1(n, &df->sk)); 629 630 PetscCall(PetscFree(df->yk)); 631 PetscCall(PetscMalloc1(n, &df->yk)); 632 633 PetscCall(PetscFree(df->ipt)); 634 PetscCall(PetscMalloc1(n, &df->ipt)); 635 636 PetscCall(PetscFree(df->ipt2)); 637 PetscCall(PetscMalloc1(n, &df->ipt2)); 638 639 PetscCall(PetscFree(df->uv)); 640 PetscCall(PetscMalloc1(n, &df->uv)); 641 PetscFunctionReturn(PETSC_SUCCESS); 642 } 643 644 static PetscErrorCode destroy_df_solver(TAO_DF *df) 645 { 646 PetscInt i; 647 648 PetscFunctionBegin; 649 PetscCall(PetscFree(df->f)); 650 PetscCall(PetscFree(df->a)); 651 PetscCall(PetscFree(df->l)); 652 PetscCall(PetscFree(df->u)); 653 PetscCall(PetscFree(df->x)); 654 655 for (i = 0; i < df->cur_num_cp; i++) PetscCall(PetscFree(df->Q[i])); 656 PetscCall(PetscFree(df->Q)); 657 PetscCall(PetscFree(df->ipt)); 658 PetscCall(PetscFree(df->ipt2)); 659 PetscCall(PetscFree(df->uv)); 660 PetscCall(PetscFree(df->g)); 661 PetscCall(PetscFree(df->y)); 662 PetscCall(PetscFree(df->tempv)); 663 PetscCall(PetscFree(df->d)); 664 PetscCall(PetscFree(df->Qd)); 665 PetscCall(PetscFree(df->t)); 666 PetscCall(PetscFree(df->xplus)); 667 PetscCall(PetscFree(df->tplus)); 668 PetscCall(PetscFree(df->sk)); 669 PetscCall(PetscFree(df->yk)); 670 PetscFunctionReturn(PETSC_SUCCESS); 671 } 672 673 /* Piecewise linear monotone target function for the Dai-Fletcher projector */ 674 static PetscReal phi(PetscReal *x, PetscInt n, PetscReal lambda, PetscReal *a, PetscReal b, PetscReal *c, PetscReal *l, PetscReal *u) 675 { 676 PetscReal r = 0.0; 677 PetscInt i; 678 679 for (i = 0; i < n; i++) { 680 x[i] = -c[i] + lambda * a[i]; 681 if (x[i] > u[i]) x[i] = u[i]; 682 else if (x[i] < l[i]) x[i] = l[i]; 683 r += a[i] * x[i]; 684 } 685 return r - b; 686 } 687 688 /** Modified Dai-Fletcher QP projector solves the problem: 689 * 690 * minimise 0.5*x'*x - c'*x 691 * subj to a'*x = b 692 * l \leq x \leq u 693 * 694 * \param c The point to be projected onto feasible set 695 */ 696 static PetscInt project(PetscInt n, PetscReal *a, PetscReal b, PetscReal *c, PetscReal *l, PetscReal *u, PetscReal *x, PetscReal *lam_ext, TAO_DF *df) 697 { 698 PetscReal lambda, lambdal, lambdau, dlambda, lambda_new; 699 PetscReal r, rl, ru, s; 700 PetscInt innerIter; 701 PetscBool nonNegativeSlack = PETSC_FALSE; 702 703 *lam_ext = 0; 704 lambda = 0; 705 dlambda = 0.5; 706 innerIter = 1; 707 708 /* \phi(x;lambda) := 0.5*x'*x + c'*x - lambda*(a'*x-b) 709 * 710 * Optimality conditions for \phi: 711 * 712 * 1. lambda <= 0 713 * 2. r <= 0 714 * 3. r*lambda == 0 715 */ 716 717 /* Bracketing Phase */ 718 r = phi(x, n, lambda, a, b, c, l, u); 719 720 if (nonNegativeSlack) { 721 /* inequality constraint, i.e., with \xi >= 0 constraint */ 722 if (r < TOL_R) return PETSC_SUCCESS; 723 } else { 724 /* equality constraint ,i.e., without \xi >= 0 constraint */ 725 if (PetscAbsReal(r) < TOL_R) return PETSC_SUCCESS; 726 } 727 728 if (r < 0.0) { 729 lambdal = lambda; 730 rl = r; 731 lambda = lambda + dlambda; 732 r = phi(x, n, lambda, a, b, c, l, u); 733 while (r < 0.0 && dlambda < BMRM_INFTY) { 734 lambdal = lambda; 735 s = rl / r - 1.0; 736 if (s < 0.1) s = 0.1; 737 dlambda = dlambda + dlambda / s; 738 lambda = lambda + dlambda; 739 rl = r; 740 r = phi(x, n, lambda, a, b, c, l, u); 741 } 742 lambdau = lambda; 743 ru = r; 744 } else { 745 lambdau = lambda; 746 ru = r; 747 lambda = lambda - dlambda; 748 r = phi(x, n, lambda, a, b, c, l, u); 749 while (r > 0.0 && dlambda > -BMRM_INFTY) { 750 lambdau = lambda; 751 s = ru / r - 1.0; 752 if (s < 0.1) s = 0.1; 753 dlambda = dlambda + dlambda / s; 754 lambda = lambda - dlambda; 755 ru = r; 756 r = phi(x, n, lambda, a, b, c, l, u); 757 } 758 lambdal = lambda; 759 rl = r; 760 } 761 762 PetscCheck(PetscAbsReal(dlambda) <= BMRM_INFTY, PETSC_COMM_SELF, PETSC_ERR_PLIB, "L2N2_DaiFletcherPGM detected Infeasible QP problem!"); 763 764 if (ru == 0) return innerIter; 765 766 /* Secant Phase */ 767 s = 1.0 - rl / ru; 768 dlambda = dlambda / s; 769 lambda = lambdau - dlambda; 770 r = phi(x, n, lambda, a, b, c, l, u); 771 772 while (PetscAbsReal(r) > TOL_R && dlambda > TOL_LAM * (1.0 + PetscAbsReal(lambda)) && innerIter < df->maxProjIter) { 773 innerIter++; 774 if (r > 0.0) { 775 if (s <= 2.0) { 776 lambdau = lambda; 777 ru = r; 778 s = 1.0 - rl / ru; 779 dlambda = (lambdau - lambdal) / s; 780 lambda = lambdau - dlambda; 781 } else { 782 s = ru / r - 1.0; 783 if (s < 0.1) s = 0.1; 784 dlambda = (lambdau - lambda) / s; 785 lambda_new = 0.75 * lambdal + 0.25 * lambda; 786 if (lambda_new < (lambda - dlambda)) lambda_new = lambda - dlambda; 787 lambdau = lambda; 788 ru = r; 789 lambda = lambda_new; 790 s = (lambdau - lambdal) / (lambdau - lambda); 791 } 792 } else { 793 if (s >= 2.0) { 794 lambdal = lambda; 795 rl = r; 796 s = 1.0 - rl / ru; 797 dlambda = (lambdau - lambdal) / s; 798 lambda = lambdau - dlambda; 799 } else { 800 s = rl / r - 1.0; 801 if (s < 0.1) s = 0.1; 802 dlambda = (lambda - lambdal) / s; 803 lambda_new = 0.75 * lambdau + 0.25 * lambda; 804 if (lambda_new > (lambda + dlambda)) lambda_new = lambda + dlambda; 805 lambdal = lambda; 806 rl = r; 807 lambda = lambda_new; 808 s = (lambdau - lambdal) / (lambdau - lambda); 809 } 810 } 811 r = phi(x, n, lambda, a, b, c, l, u); 812 } 813 814 *lam_ext = lambda; 815 if (innerIter >= df->maxProjIter) PetscCall(PetscInfo(NULL, "WARNING: DaiFletcher max iterations\n")); 816 return innerIter; 817 } 818