Lines Matching defs:cg
15 TAO_BNCG *cg = (TAO_BNCG *)tao->data;
18 PetscCall(ISDestroy(&cg->inactive_old));
19 if (cg->inactive_idx) {
20 PetscCall(ISDuplicate(cg->inactive_idx, &cg->inactive_old));
21 PetscCall(ISCopy(cg->inactive_idx, cg->inactive_old));
25 PetscCall(ISDestroy(&cg->inactive_idx));
26 PetscCall(VecWhichInactive(tao->XL, tao->solution, cg->unprojected_gradient, tao->XU, PETSC_TRUE, &cg->inactive_idx));
27 PetscCall(ISDestroy(&cg->active_idx));
28 PetscCall(ISComplementVec(cg->inactive_idx, tao->solution, &cg->active_idx));
32 PetscCall(VecCopy(cg->unprojected_gradient, cg->W));
33 PetscCall(VecScale(cg->W, -1.0));
34 PetscCall(TaoEstimateActiveBounds(tao->solution, tao->XL, tao->XU, cg->unprojected_gradient, cg->W, cg->work, cg->as_step, &cg->as_tol, &cg->active_lower, &cg->active_upper, &cg->active_fixed, &cg->active_idx, &cg->inactive_idx));
44 TAO_BNCG *cg = (TAO_BNCG *)tao->data;
49 if (cg->active_idx) PetscCall(VecISSet(step, cg->active_idx, 0.0));
52 PetscCall(TaoBoundStep(tao->solution, tao->XL, tao->XU, cg->active_lower, cg->active_upper, cg->active_fixed, 1.0, step));
62 TAO_BNCG *cg = (TAO_BNCG *)tao->data;
74 if (nDiff > 0 || !tao->recycle) PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &cg->f, cg->unprojected_gradient));
75 PetscCall(VecNorm(cg->unprojected_gradient, NORM_2, &gnorm));
76 PetscCheck(!PetscIsInfOrNanReal(cg->f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated infinity or NaN");
79 PetscCall(TaoBNCGEstimateActiveSet(tao, cg->as_type));
82 PetscCall(VecCopy(cg->unprojected_gradient, tao->gradient));
83 if (cg->active_idx) PetscCall(VecISSet(tao->gradient, cg->active_idx, 0.0));
89 cg->ls_fails = cg->descent_error = 0;
90 cg->resets = -1;
91 cg->skipped_updates = cg->pure_gd_steps = 0;
92 cg->iter_quad = 0;
97 PetscCall(VecFischer(tao->solution, cg->unprojected_gradient, tao->XL, tao->XU, cg->W));
98 PetscCall(VecNorm(cg->W, NORM_2, &resnorm));
100 PetscCall(TaoLogConvergenceHistory(tao, cg->f, resnorm, 0.0, tao->ksp_its));
101 PetscCall(TaoMonitor(tao, tao->niter, cg->f, resnorm, 0.0, step));
114 PetscCall(TaoComputeObjective(tao, tao->solution, &cg->f));
123 TAO_BNCG *cg = (TAO_BNCG *)tao->data;
128 if (!cg->W) PetscCall(VecDuplicate(tao->solution, &cg->W));
129 if (!cg->work) PetscCall(VecDuplicate(tao->solution, &cg->work));
130 if (!cg->sk) PetscCall(VecDuplicate(tao->solution, &cg->sk));
131 if (!cg->yk) PetscCall(VecDuplicate(tao->gradient, &cg->yk));
132 if (!cg->X_old) PetscCall(VecDuplicate(tao->solution, &cg->X_old));
133 if (!cg->G_old) PetscCall(VecDuplicate(tao->gradient, &cg->G_old));
134 if (cg->diag_scaling) {
135 PetscCall(VecDuplicate(tao->solution, &cg->d_work));
136 PetscCall(VecDuplicate(tao->solution, &cg->y_work));
137 PetscCall(VecDuplicate(tao->solution, &cg->g_work));
139 if (!cg->unprojected_gradient) PetscCall(VecDuplicate(tao->gradient, &cg->unprojected_gradient));
140 if (!cg->unprojected_gradient_old) PetscCall(VecDuplicate(tao->gradient, &cg->unprojected_gradient_old));
141 PetscCall(MatLMVMAllocate(cg->B, cg->sk, cg->yk));
142 if (cg->pc) PetscCall(MatLMVMSetJ0(cg->B, cg->pc));
148 TAO_BNCG *cg = (TAO_BNCG *)tao->data;
152 PetscCall(VecDestroy(&cg->W));
153 PetscCall(VecDestroy(&cg->work));
154 PetscCall(VecDestroy(&cg->X_old));
155 PetscCall(VecDestroy(&cg->G_old));
156 PetscCall(VecDestroy(&cg->unprojected_gradient));
157 PetscCall(VecDestroy(&cg->unprojected_gradient_old));
158 PetscCall(VecDestroy(&cg->g_work));
159 PetscCall(VecDestroy(&cg->d_work));
160 PetscCall(VecDestroy(&cg->y_work));
161 PetscCall(VecDestroy(&cg->sk));
162 PetscCall(VecDestroy(&cg->yk));
164 PetscCall(ISDestroy(&cg->active_lower));
165 PetscCall(ISDestroy(&cg->active_upper));
166 PetscCall(ISDestroy(&cg->active_fixed));
167 PetscCall(ISDestroy(&cg->active_idx));
168 PetscCall(ISDestroy(&cg->inactive_idx));
169 PetscCall(ISDestroy(&cg->inactive_old));
170 PetscCall(ISDestroy(&cg->new_inactives));
171 PetscCall(MatDestroy(&cg->B));
172 if (cg->pc) PetscCall(MatDestroy(&cg->pc));
179 TAO_BNCG *cg = (TAO_BNCG *)tao->data;
183 PetscCall(PetscOptionsEnum("-tao_bncg_type", "CG update formula", "TaoBNCGTypes", TaoBNCGTypes, (PetscEnum)cg->cg_type, (PetscEnum *)&cg->cg_type, NULL));
184 if (cg->cg_type != TAO_BNCG_SSML_BFGS) cg->alpha = -1.0; /* Setting defaults for non-BFGS methods. User can change it below. */
185 if (TAO_BNCG_GD == cg->cg_type) {
186 cg->cg_type = TAO_BNCG_PCGD;
188 cg->unscaled_restart = PETSC_TRUE;
189 cg->diag_scaling = PETSC_FALSE;
191 PetscCall(PetscOptionsReal("-tao_bncg_hz_eta", "(developer) cutoff tolerance for HZ", "", cg->hz_eta, &cg->hz_eta, NULL));
192 PetscCall(PetscOptionsReal("-tao_bncg_eps", "(developer) cutoff value for restarts", "", cg->epsilon, &cg->epsilon, NULL));
193 PetscCall(PetscOptionsReal("-tao_bncg_dk_eta", "(developer) cutoff tolerance for DK", "", cg->dk_eta, &cg->dk_eta, NULL));
194 PetscCall(PetscOptionsReal("-tao_bncg_xi", "(developer) Parameter in the KD method", "", cg->xi, &cg->xi, NULL));
195 PetscCall(PetscOptionsReal("-tao_bncg_theta", "(developer) update parameter for the Broyden method", "", cg->theta, &cg->theta, NULL));
196 PetscCall(PetscOptionsReal("-tao_bncg_hz_theta", "(developer) parameter for the HZ (2006) method", "", cg->hz_theta, &cg->hz_theta, NULL));
197 PetscCall(PetscOptionsReal("-tao_bncg_alpha", "(developer) parameter for the scalar scaling", "", cg->alpha, &cg->alpha, NULL));
198 PetscCall(PetscOptionsReal("-tao_bncg_bfgs_scale", "(developer) update parameter for bfgs/brdn CG methods", "", cg->bfgs_scale, &cg->bfgs_scale, NULL));
199 PetscCall(PetscOptionsReal("-tao_bncg_dfp_scale", "(developer) update parameter for bfgs/brdn CG methods", "", cg->dfp_scale, &cg->dfp_scale, NULL));
200 PetscCall(PetscOptionsBool("-tao_bncg_diag_scaling", "Enable diagonal Broyden-like preconditioning", "", cg->diag_scaling, &cg->diag_scaling, NULL));
201 PetscCall(PetscOptionsBool("-tao_bncg_dynamic_restart", "(developer) use dynamic restarts as in HZ, DK, KD", "", cg->use_dynamic_restart, &cg->use_dynamic_restart, NULL));
202 PetscCall(PetscOptionsBool("-tao_bncg_unscaled_restart", "(developer) use unscaled gradient restarts", "", cg->unscaled_restart, &cg->unscaled_restart, NULL));
203 PetscCall(PetscOptionsReal("-tao_bncg_zeta", "(developer) Free parameter for the Kou-Dai method", "", cg->zeta, &cg->zeta, NULL));
204 PetscCall(PetscOptionsInt("-tao_bncg_min_quad", "(developer) Number of iterations with approximate quadratic behavior needed for restart", "", cg->min_quad, &cg->min_quad, NULL));
205 PetscCall(PetscOptionsInt("-tao_bncg_min_restart_num", "(developer) Number of iterations between restarts (times dimension)", "", cg->min_restart_num, &cg->min_restart_num, NULL));
206 PetscCall(PetscOptionsBool("-tao_bncg_spaced_restart", "(developer) Enable regular steepest descent restarting every fixed number of iterations", "", cg->spaced_restart, &cg->spaced_restart, NULL));
207 PetscCall(PetscOptionsBool("-tao_bncg_no_scaling", "Disable all scaling except in restarts", "", cg->no_scaling, &cg->no_scaling, NULL));
208 if (cg->no_scaling) {
209 cg->diag_scaling = PETSC_FALSE;
210 cg->alpha = -1.0;
212 if (cg->alpha == -1.0 && cg->cg_type == TAO_BNCG_KD && !cg->diag_scaling) { /* Some more default options that appear to be good. */
213 cg->neg_xi = PETSC_TRUE;
215 PetscCall(PetscOptionsBool("-tao_bncg_neg_xi", "(developer) Use negative xi when it might be a smaller descent direction than necessary", "", cg->neg_xi, &cg->neg_xi, NULL));
216 PetscCall(PetscOptionsEList("-tao_bncg_as_type", "active set estimation method", "", CG_AS_TYPE, CG_AS_SIZE, CG_AS_TYPE[cg->as_type], &cg->as_type, NULL));
217 PetscCall(PetscOptionsReal("-tao_bncg_as_tol", "(developer) initial tolerance used when estimating actively bounded variables", "", cg->as_tol, &cg->as_tol, NULL));
218 PetscCall(PetscOptionsReal("-tao_bncg_as_step", "(developer) step length used when estimating actively bounded variables", "", cg->as_step, &cg->as_step, NULL));
219 PetscCall(PetscOptionsReal("-tao_bncg_delta_min", "(developer) minimum scaling factor used for scaled gradient restarts", "", cg->delta_min, &cg->delta_min, NULL));
220 PetscCall(PetscOptionsReal("-tao_bncg_delta_max", "(developer) maximum scaling factor used for scaled gradient restarts", "", cg->delta_max, &cg->delta_max, NULL));
223 PetscCall(MatSetOptionsPrefix(cg->B, ((PetscObject)tao)->prefix));
224 PetscCall(MatAppendOptionsPrefix(cg->B, "tao_bncg_"));
225 PetscCall(MatSetFromOptions(cg->B));
232 TAO_BNCG *cg = (TAO_BNCG *)tao->data;
238 PetscCall(PetscViewerASCIIPrintf(viewer, "CG Type: %s\n", TaoBNCGTypes[cg->cg_type]));
239 PetscCall(PetscViewerASCIIPrintf(viewer, "Skipped Stepdirection Updates: %" PetscInt_FMT "\n", cg->skipped_updates));
240 PetscCall(PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %" PetscInt_FMT "\n", cg->resets));
241 PetscCall(PetscViewerASCIIPrintf(viewer, "Pure gradient steps: %" PetscInt_FMT "\n", cg->pure_gd_steps));
242 PetscCall(PetscViewerASCIIPrintf(viewer, "Not a descent direction: %" PetscInt_FMT "\n", cg->descent_error));
243 PetscCall(PetscViewerASCIIPrintf(viewer, "Line search fails: %" PetscInt_FMT "\n", cg->ls_fails));
244 if (cg->diag_scaling) {
248 PetscCall(MatView(cg->B, viewer));
289 . -tao_bncg_type <taocg_type> - cg formula
333 TAO_BNCG *cg;
357 PetscCall(PetscNew(&cg));
358 tao->data = (void *)cg;
360 PetscCall(MatCreate(PetscObjectComm((PetscObject)tao), &cg->B));
361 PetscCall(PetscObjectIncrementTabLevel((PetscObject)cg->B, (PetscObject)tao, 1));
362 PetscCall(MatSetType(cg->B, MATLMVMDIAGBROYDEN));
364 cg->pc = NULL;
366 cg->dk_eta = 0.5;
367 cg->hz_eta = 0.4;
368 cg->dynamic_restart = PETSC_FALSE;
369 cg->unscaled_restart = PETSC_FALSE;
370 cg->no_scaling = PETSC_FALSE;
371 cg->delta_min = 1e-7;
372 cg->delta_max = 100;
373 cg->theta = 1.0;
374 cg->hz_theta = 1.0;
375 cg->dfp_scale = 1.0;
376 cg->bfgs_scale = 1.0;
377 cg->zeta = 0.1;
378 cg->min_quad = 6;
379 cg->min_restart_num = 6; /* As in CG_DESCENT and KD2015*/
380 cg->xi = 1.0;
381 cg->neg_xi = PETSC_TRUE;
382 cg->spaced_restart = PETSC_FALSE;
383 cg->tol_quad = 1e-8;
384 cg->as_step = 0.001;
385 cg->as_tol = 0.001;
386 cg->eps_23 = PetscPowReal(PETSC_MACHINE_EPSILON, 2.0 / 3.0); /* Just a little tighter*/
387 cg->as_type = CG_AS_BERTSEKAS;
388 cg->cg_type = TAO_BNCG_SSML_BFGS;
389 cg->alpha = 1.0;
390 cg->diag_scaling = PETSC_TRUE;
396 TAO_BNCG *cg = (TAO_BNCG *)tao->data;
400 ++cg->resets;
401 scaling = 2.0 * PetscMax(1.0, PetscAbsScalar(cg->f)) / PetscMax(gnormsq, cg->eps_23);
402 scaling = PetscMin(cg->delta_max, PetscMax(cg->delta_min, scaling));
403 if (cg->unscaled_restart) {
405 ++cg->pure_gd_steps;
409 if (cg->diag_scaling) PetscCall(MatLMVMReset(cg->B, PETSC_FALSE));
415 TAO_BNCG *cg = (TAO_BNCG *)tao->data;
419 if (cg->f < cg->min_quad / 10) {
423 quadinterp = 2.0 * (cg->f - fold) / (stepsize * (gd + gd_old));
424 if (PetscAbs(quadinterp - 1.0) < cg->tol_quad) ++cg->iter_quad;
426 cg->iter_quad = 0;
429 if (cg->iter_quad >= cg->min_quad) {
430 cg->iter_quad = 0;
438 TAO_BNCG *cg = (TAO_BNCG *)tao->data;
448 PetscCall(VecWAXPY(cg->yk, -1.0, cg->G_old, tao->gradient));
449 PetscCall(VecNorm(cg->yk, NORM_2, &ynorm));
451 PetscCall(VecDot(cg->yk, tao->stepdirection, &dk_yk));
454 ++cg->skipped_updates;
456 if (cg->spaced_restart) {
458 if (tao->niter % (dim * cg->min_restart_num)) cg_restart = PETSC_TRUE;
462 if (cg->spaced_restart) {
467 if (cg->diag_scaling) PetscCall(MatLMVMUpdate(cg->B, tao->solution, tao->gradient));
495 PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
496 PetscCall(VecAXPBY(tao->stepdirection, -1.0, 0.0, cg->g_work));
498 switch (cg->cg_type) {
500 if (!cg->diag_scaling) {
501 if (!cg->no_scaling) {
502 cg->sts = step * step * dnorm * dnorm;
503 PetscCall(TaoBNCGComputeScalarScaling(ynorm2, step * dk_yk, cg->sts, &tau_k, cg->alpha));
506 ++cg->pure_gd_steps;
510 PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
511 PetscCall(VecAXPBY(tao->stepdirection, -1.0, 0.0, cg->g_work));
517 if (!cg->diag_scaling) {
518 cg->sts = step * step * dnorm * dnorm;
519 PetscCall(VecDot(cg->yk, tao->gradient, &gkp1_yk));
520 PetscCall(TaoBNCGComputeScalarScaling(ynorm2, step * dk_yk, cg->sts, &tau_k, cg->alpha));
524 PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
525 PetscCall(VecDot(cg->yk, cg->g_work, &gkp1_yk));
527 PetscCall(VecAXPBY(tao->stepdirection, -1.0, beta, cg->g_work));
532 PetscCall(VecDot(cg->G_old, cg->G_old, &gnorm2_old));
533 PetscCall(VecWAXPY(cg->yk, -1.0, cg->G_old, tao->gradient));
534 PetscCall(VecNorm(cg->yk, NORM_2, &ynorm));
536 PetscCall(VecDot(cg->yk, tao->stepdirection, &dk_yk));
537 if (!cg->diag_scaling) {
538 PetscCall(TaoBNCGComputeScalarScaling(ynorm2, step * dk_yk, step * step * dnorm * dnorm, &tau_k, cg->alpha));
542 PetscCall(VecDot(cg->G_old, cg->g_work, &gnorm2_old)); /* Before it's updated */
543 PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
544 PetscCall(VecDot(tao->gradient, cg->g_work, &tmp));
546 PetscCall(VecAXPBY(tao->stepdirection, -1.0, beta, cg->g_work));
552 if (!cg->diag_scaling) {
553 PetscCall(VecDot(cg->G_old, cg->G_old, &gnorm2_old));
554 PetscCall(VecDot(cg->yk, tao->gradient, &gkp1_yk));
555 PetscCall(TaoBNCGComputeScalarScaling(ynorm2, step * dk_yk, snorm * snorm, &tau_k, cg->alpha));
559 PetscCall(VecDot(cg->G_old, cg->g_work, &gnorm2_old));
560 PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
561 PetscCall(VecDot(cg->g_work, cg->yk, &gkp1_yk));
563 PetscCall(VecAXPBY(tao->stepdirection, -1.0, beta, cg->g_work));
568 PetscCall(VecWAXPY(cg->yk, -1.0, cg->G_old, tao->gradient));
569 PetscCall(VecNorm(cg->yk, NORM_2, &ynorm));
571 if (!cg->diag_scaling) {
572 PetscCall(VecDot(cg->G_old, cg->G_old, &gnorm2_old));
573 PetscCall(VecDot(cg->yk, tao->gradient, &gkp1_yk));
574 PetscCall(TaoBNCGComputeScalarScaling(ynorm2, step * dk_yk, snorm * snorm, &tau_k, cg->alpha));
579 PetscCall(VecDot(cg->G_old, cg->g_work, &gnorm2_old)); /* Old gtDg */
580 PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
581 PetscCall(VecDot(cg->g_work, cg->yk, &gkp1_yk));
584 PetscCall(VecAXPBY(tao->stepdirection, -1.0, beta, cg->g_work));
591 if (!cg->diag_scaling) {
593 PetscCall(VecDot(cg->G_old, tao->stepdirection, &gd_old));
594 PetscCall(TaoBNCGComputeScalarScaling(ynorm2, step * dk_yk, cg->yts, &tau_k, cg->alpha));
598 PetscCall(MatMult(cg->B, tao->stepdirection, cg->d_work));
599 PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
600 PetscCall(VecDot(cg->g_work, tao->gradient, >Dg));
601 PetscCall(VecDot(tao->stepdirection, cg->G_old, &gd_old));
602 PetscCall(VecDot(cg->d_work, cg->g_work, &dk_yk));
605 PetscCall(VecScale(cg->d_work, beta));
606 PetscCall(VecWAXPY(tao->stepdirection, -1.0, cg->g_work, cg->d_work));
614 PetscCall(VecDot(cg->G_old, tao->stepdirection, &gd_old));
615 PetscCall(VecWAXPY(cg->sk, -1.0, cg->X_old, tao->solution));
617 cg->yts = step * dk_yk;
618 if (cg->use_dynamic_restart) PetscCall(TaoBNCGCheckDynamicRestart(tao, step, gd, gd_old, &cg->dynamic_restart, fold));
619 if (cg->dynamic_restart) {
622 if (!cg->diag_scaling) {
623 PetscCall(VecDot(cg->yk, tao->gradient, &gkp1_yk));
624 PetscCall(TaoBNCGComputeScalarScaling(ynorm2, cg->yts, snorm * snorm, &tau_k, cg->alpha));
625 /* Supplying cg->alpha = -1.0 will give the CG_DESCENT 5.3 special case of tau_k = 1.0 */
629 beta = PetscMax(PetscMax(beta, cg->hz_eta * tau_k * gd_old / (dnorm * dnorm)), cg->dk_eta * tau_k * gd / (dnorm * dnorm));
634 cg->yty = ynorm2;
635 cg->sts = snorm * snorm;
637 PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
638 PetscCall(MatSolve(cg->B, cg->yk, cg->y_work));
639 PetscCall(MatSolve(cg->B, tao->stepdirection, cg->d_work));
641 PetscCall(VecDot(cg->yk, cg->g_work, &gkp1_yk));
644 PetscCall(VecDot(cg->yk, cg->y_work, &tau_k));
647 beta = gkp1_yk / dk_yk + cg->hz_theta * tau_k; /* HZ; (1.15) from DK 2013 */
649 PetscCall(VecDot(cg->G_old, cg->d_work, &gd_old));
650 PetscCall(VecDot(tao->stepdirection, cg->g_work, &gd));
651 beta = PetscMax(PetscMax(beta, cg->hz_eta * gd_old / (dnorm * dnorm)), cg->dk_eta * gd / (dnorm * dnorm));
653 PetscCall(VecAXPBY(tao->stepdirection, -1.0, beta, cg->g_work));
662 PetscCall(VecDot(cg->G_old, tao->stepdirection, &gd_old));
663 PetscCall(VecWAXPY(cg->sk, -1.0, cg->X_old, tao->solution));
665 cg->yts = dk_yk * step;
666 if (!cg->diag_scaling) {
667 PetscCall(VecDot(cg->yk, tao->gradient, &gkp1_yk));
668 PetscCall(TaoBNCGComputeScalarScaling(ynorm2, cg->yts, snorm * snorm, &tau_k, cg->alpha));
669 /* Use cg->alpha = -1.0 to get tau_k = 1.0 as in CG_DESCENT 5.3 */
672 beta = PetscMax(PetscMax(beta, cg->hz_eta * tau_k * gd_old / (dnorm * dnorm)), cg->dk_eta * tau_k * gd / (dnorm * dnorm));
674 PetscCall(VecAXPBYPCZ(tao->stepdirection, -tau_k, 0.0, beta, tao->gradient, cg->yk));
677 cg->yty = ynorm2;
678 cg->sts = snorm * snorm;
679 PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
680 PetscCall(MatSolve(cg->B, cg->yk, cg->y_work));
681 PetscCall(MatSolve(cg->B, tao->stepdirection, cg->d_work));
683 PetscCall(VecDot(cg->yk, cg->g_work, &gkp1_yk));
684 PetscCall(VecDot(cg->yk, cg->y_work, &tau_k));
690 PetscCall(VecDot(cg->y_work, tao->stepdirection, &dk_yk));
692 PetscCall(VecDot(tao->stepdirection, cg->g_work, &gd));
693 PetscCall(VecDot(cg->G_old, cg->d_work, &gd_old));
694 beta = PetscMax(PetscMax(beta, cg->hz_eta * gd_old / (dnorm * dnorm)), cg->dk_eta * gd / (dnorm * dnorm));
696 PetscCall(VecAXPBY(tao->stepdirection, -1.0, beta, cg->g_work));
704 PetscCall(VecDot(cg->G_old, tao->stepdirection, &gd_old));
705 PetscCall(VecWAXPY(cg->sk, -1.0, cg->X_old, tao->solution));
707 cg->yts = dk_yk * step;
708 if (cg->use_dynamic_restart) PetscCall(TaoBNCGCheckDynamicRestart(tao, step, gd, gd_old, &cg->dynamic_restart, fold));
709 if (cg->dynamic_restart) {
712 if (!cg->diag_scaling) {
713 PetscCall(VecDot(cg->yk, tao->gradient, &gkp1_yk));
714 PetscCall(TaoBNCGComputeScalarScaling(ynorm2, cg->yts, snorm * snorm, &tau_k, cg->alpha));
716 if (beta < cg->zeta * tau_k * gd / (dnorm * dnorm)) /* 0.1 is KD's zeta parameter */
718 beta = cg->zeta * tau_k * gd / (dnorm * dnorm);
721 if (gkp1_yk < 0 && cg->neg_xi) gamma = -1.0 * gd / dk_yk;
724 else gamma = cg->xi * gd / dk_yk;
727 PetscCall(VecAXPBYPCZ(tao->stepdirection, -tau_k, gamma * tau_k, beta, tao->gradient, cg->yk));
730 cg->yty = ynorm2;
731 cg->sts = snorm * snorm;
732 PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
733 PetscCall(MatSolve(cg->B, cg->yk, cg->y_work));
735 PetscCall(VecDot(cg->yk, cg->g_work, &gkp1D_yk));
739 PetscCall(VecDot(cg->yk, cg->y_work, &tau_k));
744 PetscCall(VecDot(tao->stepdirection, cg->g_work, &tmp));
745 if (cg->neg_xi) {
748 else gamma = cg->xi * gd / dk_yk;
749 if (beta < cg->zeta * tmp / (dnorm * dnorm)) {
750 beta = cg->zeta * tmp / (dnorm * dnorm);
754 if (beta < cg->zeta * tmp / (dnorm * dnorm)) {
755 beta = cg->zeta * tmp / (dnorm * dnorm);
757 } else gamma = cg->xi * gd / dk_yk;
760 PetscCall(VecAXPBY(tao->stepdirection, -1.0, beta, cg->g_work));
761 PetscCall(VecAXPY(tao->stepdirection, gamma, cg->y_work));
770 PetscCall(VecWAXPY(cg->sk, -1.0, cg->X_old, tao->solution));
772 cg->yts = dk_yk * step;
773 cg->yty = ynorm2;
774 cg->sts = snorm * snorm;
775 if (!cg->diag_scaling) {
776 PetscCall(VecDot(cg->yk, tao->gradient, &gkp1_yk));
777 PetscCall(TaoBNCGComputeScalarScaling(cg->yty, cg->yts, cg->sts, &tau_k, cg->alpha));
779 beta = tau_k * (gkp1_yk / dk_yk - cg->yty * gd / (dk_yk * dk_yk)) - step * tmp;
781 PetscCall(VecAXPBYPCZ(tao->stepdirection, -tau_k, tmp * tau_k, beta, tao->gradient, cg->yk));
784 PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
785 PetscCall(MatSolve(cg->B, cg->yk, cg->y_work));
787 PetscCall(VecDot(cg->g_work, cg->yk, &gkp1_yk));
788 PetscCall(VecDot(cg->y_work, cg->yk, &tmp));
793 PetscCall(VecAXPBYPCZ(tao->stepdirection, -1.0, gamma, beta, cg->g_work, cg->y_work));
799 PetscCall(VecWAXPY(cg->sk, -1.0, cg->X_old, tao->solution));
801 cg->yts = dk_yk * step;
802 cg->yty = ynorm2;
803 cg->sts = snorm * snorm;
804 if (!cg->diag_scaling) {
806 PetscCall(TaoBNCGComputeScalarScaling(cg->yty, cg->yts, cg->sts, &tau_k, cg->alpha));
807 PetscCall(VecDot(cg->yk, tao->gradient, &gkp1_yk));
808 tau_k = cg->dfp_scale * tau_k;
809 tmp = tau_k * gkp1_yk / cg->yty;
812 PetscCall(VecAXPBYPCZ(tao->stepdirection, -tau_k, tmp, beta, tao->gradient, cg->yk));
815 PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
816 PetscCall(MatSolve(cg->B, cg->yk, cg->y_work));
818 PetscCall(VecDot(cg->g_work, cg->yk, &gkp1_yk));
819 PetscCall(VecDot(cg->y_work, cg->yk, &tmp));
824 PetscCall(VecAXPBYPCZ(tao->stepdirection, -1.0, gamma, beta, cg->g_work, cg->y_work));
830 PetscCall(VecWAXPY(cg->sk, -1.0, cg->X_old, tao->solution));
832 cg->yts = step * dk_yk;
833 cg->yty = ynorm2;
834 cg->sts = snorm * snorm;
835 if (!cg->diag_scaling) {
837 PetscCall(TaoBNCGComputeScalarScaling(cg->yty, step * dk_yk, snorm * snorm, &tau_bfgs, cg->bfgs_scale));
838 PetscCall(TaoBNCGComputeScalarScaling(cg->yty, step * dk_yk, snorm * snorm, &tau_dfp, cg->dfp_scale));
839 PetscCall(VecDot(cg->yk, tao->gradient, &gkp1_yk));
840 tau_k = cg->theta * tau_bfgs + (1.0 - cg->theta) * tau_dfp;
843 tmp = cg->theta * tau_bfgs * gd / dk_yk + (1 - cg->theta) * tau_dfp * gkp1_yk / cg->yty;
844 beta = cg->theta * tau_bfgs * (gkp1_yk / dk_yk - cg->yty * gd / (dk_yk * dk_yk)) - step * gd / dk_yk;
846 PetscCall(VecAXPBYPCZ(tao->stepdirection, -tau_k, tmp, beta, tao->gradient, cg->yk));
849 PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
850 PetscCall(MatSolve(cg->B, cg->yk, cg->y_work));
852 PetscCall(VecDot(cg->g_work, cg->yk, &gkp1_yk));
853 PetscCall(VecDot(cg->y_work, cg->yk, &tmp));
854 gamma = cg->theta * gd / dk_yk + (1 - cg->theta) * (gkp1_yk / tmp);
856 beta = cg->theta * (gkp1_yk / dk_yk - gd * tmp / (dk_yk * dk_yk)) - step * gd / dk_yk;
858 PetscCall(VecAXPBYPCZ(tao->stepdirection, -1.0, gamma, beta, cg->g_work, cg->y_work));
871 TAO_BNCG *cg = (TAO_BNCG *)tao->data;
880 PetscCall(VecCopy(tao->solution, cg->X_old));
881 PetscCall(VecCopy(tao->gradient, cg->G_old));
882 PetscCall(VecCopy(cg->unprojected_gradient, cg->unprojected_gradient_old));
886 f_old = cg->f;
892 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &cg->f, cg->unprojected_gradient, tao->stepdirection, &step, &ls_status));
897 ++cg->ls_fails;
898 if (cg->cg_type == TAO_BNCG_GD) {
904 PetscCall(VecCopy(cg->X_old, tao->solution));
905 PetscCall(VecCopy(cg->G_old, tao->gradient));
906 PetscCall(VecCopy(cg->unprojected_gradient_old, cg->unprojected_gradient));
909 cg->f = f_old;
912 if (cg->cg_type != TAO_BNCG_PCGD && cg->diag_scaling) {
917 PetscCall(TaoBNCGBoundStep(tao, cg->as_type, tao->stepdirection));
920 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &cg->f, cg->unprojected_gradient, tao->stepdirection, &step, &ls_status));
926 ++cg->ls_fails;
932 ++cg->ls_fails;
934 PetscCall(TaoBNCGBoundStep(tao, cg->as_type, tao->stepdirection));
936 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &cg->f, cg->unprojected_gradient, tao->stepdirection, &step, &ls_status));
942 ++cg->ls_fails;
955 PetscCall(VecFischer(tao->solution, cg->unprojected_gradient, tao->XL, tao->XU, cg->W));
956 PetscCall(VecNorm(cg->W, NORM_2, &resnorm));
958 PetscCall(TaoLogConvergenceHistory(tao, cg->f, resnorm, 0.0, tao->ksp_its));
959 PetscCall(TaoMonitor(tao, tao->niter, cg->f, resnorm, 0.0, step));
966 PetscCall(TaoBNCGEstimateActiveSet(tao, cg->as_type));
968 PetscCall(VecCopy(cg->unprojected_gradient, tao->gradient));
969 if (cg->active_idx) PetscCall(VecISSet(tao->gradient, cg->active_idx, 0.0));
978 PetscCall(TaoBNCGBoundStep(tao, cg->as_type, tao->stepdirection));
980 if (cg->cg_type != TAO_BNCG_GD) {
982 PetscCall(ISDestroy(&cg->new_inactives));
983 if (cg->inactive_idx && cg->inactive_old) PetscCall(ISDifference(cg->inactive_idx, cg->inactive_old, &cg->new_inactives));
985 if (cg->new_inactives) {
986 PetscCall(VecGetSubVector(tao->stepdirection, cg->new_inactives, &cg->inactive_step));
987 PetscCall(VecGetSubVector(cg->unprojected_gradient, cg->new_inactives, &cg->inactive_grad));
988 PetscCall(VecCopy(cg->inactive_grad, cg->inactive_step));
989 PetscCall(VecScale(cg->inactive_step, -1.0));
990 PetscCall(VecRestoreSubVector(tao->stepdirection, cg->new_inactives, &cg->inactive_step));
991 PetscCall(VecRestoreSubVector(cg->unprojected_gradient, cg->new_inactives, &cg->inactive_grad));
999 PetscCall(TaoBNCGBoundStep(tao, cg->as_type, tao->stepdirection));
1000 ++cg->descent_error;
1009 TAO_BNCG *cg = (TAO_BNCG *)tao->data;
1016 cg->pc = H0;
1036 TAO_BNCG *cg = (TAO_BNCG *)tao->data;
1042 *type = cg->cg_type;
1059 TAO_BNCG *cg = (TAO_BNCG *)tao->data;
1064 if (same) cg->cg_type = type;