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