Actual source code: bqnk.c
petsc-main 2021-04-20
1: #include <../src/tao/bound/impls/bqnk/bqnk.h>
2: #include <petscksp.h>
4: static const char *BQNK_INIT[64] = {"constant", "direction"};
5: static const char *BNK_UPDATE[64] = {"step", "reduction", "interpolation"};
6: static const char *BNK_AS[64] = {"none", "bertsekas"};
8: static PetscErrorCode TaoBQNKComputeHessian(Tao tao)
9: {
10: TAO_BNK *bnk = (TAO_BNK *)tao->data;
11: TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx;
13: PetscReal gnorm2, delta;
16: /* Alias the LMVM matrix into the TAO hessian */
17: if (tao->hessian) {
18: MatDestroy(&tao->hessian);
19: }
20: if (tao->hessian_pre) {
21: MatDestroy(&tao->hessian_pre);
22: }
23: PetscObjectReference((PetscObject)bqnk->B);
24: tao->hessian = bqnk->B;
25: PetscObjectReference((PetscObject)bqnk->B);
26: tao->hessian_pre = bqnk->B;
27: /* Update the Hessian with the latest solution */
28: if (bqnk->is_spd) {
29: gnorm2 = bnk->gnorm*bnk->gnorm;
30: if (gnorm2 == 0.0) gnorm2 = PETSC_MACHINE_EPSILON;
31: if (bnk->f == 0.0) {
32: delta = 2.0 / gnorm2;
33: } else {
34: delta = 2.0 * PetscAbsScalar(bnk->f) / gnorm2;
35: }
36: MatLMVMSymBroydenSetDelta(bqnk->B, delta);
37: }
38: MatLMVMUpdate(tao->hessian, tao->solution, bnk->unprojected_gradient);
39: MatLMVMResetShift(tao->hessian);
40: /* Prepare the reduced sub-matrices for the inactive set */
41: MatDestroy(&bnk->H_inactive);
42: if (bnk->active_idx) {
43: MatCreateSubMatrixVirtual(tao->hessian, bnk->inactive_idx, bnk->inactive_idx, &bnk->H_inactive);
44: PCLMVMSetIS(bqnk->pc, bnk->inactive_idx);
45: } else {
46: PetscObjectReference((PetscObject)tao->hessian);
47: bnk->H_inactive = tao->hessian;
48: PCLMVMClearIS(bqnk->pc);
49: }
50: MatDestroy(&bnk->Hpre_inactive);
51: PetscObjectReference((PetscObject)bnk->H_inactive);
52: bnk->Hpre_inactive = bnk->H_inactive;
53: return(0);
54: }
56: static PetscErrorCode TaoBQNKComputeStep(Tao tao, PetscBool shift, KSPConvergedReason *ksp_reason, PetscInt *step_type)
57: {
58: TAO_BNK *bnk = (TAO_BNK *)tao->data;
59: TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx;
63: TaoBNKComputeStep(tao, shift, ksp_reason, step_type);
64: if (*ksp_reason < 0) {
65: /* Krylov solver failed to converge so reset the LMVM matrix */
66: MatLMVMReset(bqnk->B, PETSC_FALSE);
67: MatLMVMUpdate(bqnk->B, tao->solution, bnk->unprojected_gradient);
68: }
69: return(0);
70: }
72: PetscErrorCode TaoSolve_BQNK(Tao tao)
73: {
74: TAO_BNK *bnk = (TAO_BNK *)tao->data;
75: TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx;
76: Mat_LMVM *lmvm = (Mat_LMVM*)bqnk->B->data;
77: Mat_LMVM *J0;
78: Mat_SymBrdn *diag_ctx;
79: PetscBool flg = PETSC_FALSE;
83: if (!tao->recycle) {
84: MatLMVMReset(bqnk->B, PETSC_FALSE);
85: lmvm->nresets = 0;
86: if (lmvm->J0) {
87: PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &flg);
88: if (flg) {
89: J0 = (Mat_LMVM*)lmvm->J0->data;
90: J0->nresets = 0;
91: }
92: }
93: flg = PETSC_FALSE;
94: PetscObjectTypeCompareAny((PetscObject)bqnk->B, &flg, MATLMVMSYMBROYDEN, MATLMVMSYMBADBROYDEN, MATLMVMBFGS, MATLMVMDFP, "");
95: if (flg) {
96: diag_ctx = (Mat_SymBrdn*)lmvm->ctx;
97: J0 = (Mat_LMVM*)diag_ctx->D->data;
98: J0->nresets = 0;
99: }
100: }
101: (*bqnk->solve)(tao);
102: return(0);
103: }
105: PetscErrorCode TaoSetUp_BQNK(Tao tao)
106: {
107: TAO_BNK *bnk = (TAO_BNK *)tao->data;
108: TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx;
110: PetscInt n, N;
111: PetscBool is_lmvm, is_sym, is_spd;
114: TaoSetUp_BNK(tao);
115: VecGetLocalSize(tao->solution,&n);
116: VecGetSize(tao->solution,&N);
117: MatSetSizes(bqnk->B, n, n, N, N);
118: MatLMVMAllocate(bqnk->B,tao->solution,bnk->unprojected_gradient);
119: PetscObjectBaseTypeCompare((PetscObject)bqnk->B, MATLMVM, &is_lmvm);
120: if (!is_lmvm) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "Matrix must be an LMVM-type");
121: MatGetOption(bqnk->B, MAT_SYMMETRIC, &is_sym);
122: if (!is_sym) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix must be symmetric");
123: MatGetOption(bqnk->B, MAT_SPD, &is_spd);
124: KSPGetPC(tao->ksp, &bqnk->pc);
125: PCSetType(bqnk->pc, PCLMVM);
126: PCLMVMSetMatLMVM(bqnk->pc, bqnk->B);
127: return(0);
128: }
130: static PetscErrorCode TaoSetFromOptions_BQNK(PetscOptionItems *PetscOptionsObject,Tao tao)
131: {
132: TAO_BNK *bnk = (TAO_BNK *)tao->data;
133: TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx;
137: PetscOptionsHead(PetscOptionsObject,"Quasi-Newton-Krylov method for bound constrained optimization");
138: PetscOptionsEList("-tao_bqnk_init_type", "radius initialization type", "", BQNK_INIT, BQNK_INIT_TYPES, BQNK_INIT[bnk->init_type], &bnk->init_type, NULL);
139: PetscOptionsEList("-tao_bqnk_update_type", "radius update type", "", BNK_UPDATE, BNK_UPDATE_TYPES, BNK_UPDATE[bnk->update_type], &bnk->update_type, NULL);
140: PetscOptionsEList("-tao_bqnk_as_type", "active set estimation method", "", BNK_AS, BNK_AS_TYPES, BNK_AS[bnk->as_type], &bnk->as_type, NULL);
141: PetscOptionsReal("-tao_bqnk_sval", "(developer) Hessian perturbation starting value", "", bnk->sval, &bnk->sval,NULL);
142: PetscOptionsReal("-tao_bqnk_imin", "(developer) minimum initial Hessian perturbation", "", bnk->imin, &bnk->imin,NULL);
143: PetscOptionsReal("-tao_bqnk_imax", "(developer) maximum initial Hessian perturbation", "", bnk->imax, &bnk->imax,NULL);
144: PetscOptionsReal("-tao_bqnk_imfac", "(developer) initial merit factor for Hessian perturbation", "", bnk->imfac, &bnk->imfac,NULL);
145: PetscOptionsReal("-tao_bqnk_pmin", "(developer) minimum Hessian perturbation", "", bnk->pmin, &bnk->pmin,NULL);
146: PetscOptionsReal("-tao_bqnk_pmax", "(developer) maximum Hessian perturbation", "", bnk->pmax, &bnk->pmax,NULL);
147: PetscOptionsReal("-tao_bqnk_pgfac", "(developer) Hessian perturbation growth factor", "", bnk->pgfac, &bnk->pgfac,NULL);
148: PetscOptionsReal("-tao_bqnk_psfac", "(developer) Hessian perturbation shrink factor", "", bnk->psfac, &bnk->psfac,NULL);
149: PetscOptionsReal("-tao_bqnk_pmgfac", "(developer) merit growth factor for Hessian perturbation", "", bnk->pmgfac, &bnk->pmgfac,NULL);
150: PetscOptionsReal("-tao_bqnk_pmsfac", "(developer) merit shrink factor for Hessian perturbation", "", bnk->pmsfac, &bnk->pmsfac,NULL);
151: PetscOptionsReal("-tao_bqnk_eta1", "(developer) threshold for rejecting step (-tao_bqnk_update_type reduction)", "", bnk->eta1, &bnk->eta1,NULL);
152: PetscOptionsReal("-tao_bqnk_eta2", "(developer) threshold for accepting marginal step (-tao_bqnk_update_type reduction)", "", bnk->eta2, &bnk->eta2,NULL);
153: PetscOptionsReal("-tao_bqnk_eta3", "(developer) threshold for accepting reasonable step (-tao_bqnk_update_type reduction)", "", bnk->eta3, &bnk->eta3,NULL);
154: PetscOptionsReal("-tao_bqnk_eta4", "(developer) threshold for accepting good step (-tao_bqnk_update_type reduction)", "", bnk->eta4, &bnk->eta4,NULL);
155: PetscOptionsReal("-tao_bqnk_alpha1", "(developer) radius reduction factor for rejected step (-tao_bqnk_update_type reduction)", "", bnk->alpha1, &bnk->alpha1,NULL);
156: PetscOptionsReal("-tao_bqnk_alpha2", "(developer) radius reduction factor for marginally accepted bad step (-tao_bqnk_update_type reduction)", "", bnk->alpha2, &bnk->alpha2,NULL);
157: PetscOptionsReal("-tao_bqnk_alpha3", "(developer) radius increase factor for reasonable accepted step (-tao_bqnk_update_type reduction)", "", bnk->alpha3, &bnk->alpha3,NULL);
158: PetscOptionsReal("-tao_bqnk_alpha4", "(developer) radius increase factor for good accepted step (-tao_bqnk_update_type reduction)", "", bnk->alpha4, &bnk->alpha4,NULL);
159: PetscOptionsReal("-tao_bqnk_alpha5", "(developer) radius increase factor for very good accepted step (-tao_bqnk_update_type reduction)", "", bnk->alpha5, &bnk->alpha5,NULL);
160: PetscOptionsReal("-tao_bqnk_nu1", "(developer) threshold for small line-search step length (-tao_bqnk_update_type step)", "", bnk->nu1, &bnk->nu1,NULL);
161: PetscOptionsReal("-tao_bqnk_nu2", "(developer) threshold for reasonable line-search step length (-tao_bqnk_update_type step)", "", bnk->nu2, &bnk->nu2,NULL);
162: PetscOptionsReal("-tao_bqnk_nu3", "(developer) threshold for large line-search step length (-tao_bqnk_update_type step)", "", bnk->nu3, &bnk->nu3,NULL);
163: PetscOptionsReal("-tao_bqnk_nu4", "(developer) threshold for very large line-search step length (-tao_bqnk_update_type step)", "", bnk->nu4, &bnk->nu4,NULL);
164: PetscOptionsReal("-tao_bqnk_omega1", "(developer) radius reduction factor for very small line-search step length (-tao_bqnk_update_type step)", "", bnk->omega1, &bnk->omega1,NULL);
165: PetscOptionsReal("-tao_bqnk_omega2", "(developer) radius reduction factor for small line-search step length (-tao_bqnk_update_type step)", "", bnk->omega2, &bnk->omega2,NULL);
166: PetscOptionsReal("-tao_bqnk_omega3", "(developer) radius factor for decent line-search step length (-tao_bqnk_update_type step)", "", bnk->omega3, &bnk->omega3,NULL);
167: PetscOptionsReal("-tao_bqnk_omega4", "(developer) radius increase factor for large line-search step length (-tao_bqnk_update_type step)", "", bnk->omega4, &bnk->omega4,NULL);
168: PetscOptionsReal("-tao_bqnk_omega5", "(developer) radius increase factor for very large line-search step length (-tao_bqnk_update_type step)", "", bnk->omega5, &bnk->omega5,NULL);
169: PetscOptionsReal("-tao_bqnk_mu1", "(developer) threshold for accepting very good step (-tao_bqnk_update_type interpolation)", "", bnk->mu1, &bnk->mu1,NULL);
170: PetscOptionsReal("-tao_bqnk_mu2", "(developer) threshold for accepting good step (-tao_bqnk_update_type interpolation)", "", bnk->mu2, &bnk->mu2,NULL);
171: PetscOptionsReal("-tao_bqnk_gamma1", "(developer) radius reduction factor for rejected very bad step (-tao_bqnk_update_type interpolation)", "", bnk->gamma1, &bnk->gamma1,NULL);
172: PetscOptionsReal("-tao_bqnk_gamma2", "(developer) radius reduction factor for rejected bad step (-tao_bqnk_update_type interpolation)", "", bnk->gamma2, &bnk->gamma2,NULL);
173: PetscOptionsReal("-tao_bqnk_gamma3", "(developer) radius increase factor for accepted good step (-tao_bqnk_update_type interpolation)", "", bnk->gamma3, &bnk->gamma3,NULL);
174: PetscOptionsReal("-tao_bqnk_gamma4", "(developer) radius increase factor for accepted very good step (-tao_bqnk_update_type interpolation)", "", bnk->gamma4, &bnk->gamma4,NULL);
175: PetscOptionsReal("-tao_bqnk_theta", "(developer) trust region interpolation factor (-tao_bqnk_update_type interpolation)", "", bnk->theta, &bnk->theta,NULL);
176: PetscOptionsReal("-tao_bqnk_min_radius", "(developer) lower bound on initial radius", "", bnk->min_radius, &bnk->min_radius,NULL);
177: PetscOptionsReal("-tao_bqnk_max_radius", "(developer) upper bound on radius", "", bnk->max_radius, &bnk->max_radius,NULL);
178: PetscOptionsReal("-tao_bqnk_epsilon", "(developer) tolerance used when computing actual and predicted reduction", "", bnk->epsilon, &bnk->epsilon,NULL);
179: PetscOptionsReal("-tao_bqnk_as_tol", "(developer) initial tolerance used when estimating actively bounded variables", "", bnk->as_tol, &bnk->as_tol,NULL);
180: PetscOptionsReal("-tao_bqnk_as_step", "(developer) step length used when estimating actively bounded variables", "", bnk->as_step, &bnk->as_step,NULL);
181: PetscOptionsInt("-tao_bqnk_max_cg_its", "number of BNCG iterations to take for each Newton step", "", bnk->max_cg_its, &bnk->max_cg_its,NULL);
182: PetscOptionsTail();
183: TaoSetFromOptions(bnk->bncg);
184: TaoLineSearchSetFromOptions(tao->linesearch);
185: KSPSetFromOptions(tao->ksp);
186: if (bnk->init_type == BNK_INIT_INTERPOLATION) bnk->init_type = BNK_INIT_DIRECTION;
187: MatSetFromOptions(bqnk->B);
188: MatGetOption(bqnk->B, MAT_SPD, &bqnk->is_spd);
189: return(0);
190: }
192: static PetscErrorCode TaoView_BQNK(Tao tao, PetscViewer viewer)
193: {
194: TAO_BNK *bnk = (TAO_BNK*)tao->data;
195: TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx;
197: PetscBool isascii;
200: TaoView_BNK(tao, viewer);
201: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
202: if (isascii) {
203: PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO);
204: MatView(bqnk->B, viewer);
205: PetscViewerPopFormat(viewer);
206: }
207: return(0);
208: }
210: static PetscErrorCode TaoDestroy_BQNK(Tao tao)
211: {
212: TAO_BNK *bnk = (TAO_BNK*)tao->data;
213: TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx;
217: MatDestroy(&bnk->Hpre_inactive);
218: MatDestroy(&bnk->H_inactive);
219: MatDestroy(&bqnk->B);
220: PetscFree(bnk->ctx);
221: TaoDestroy_BNK(tao);
222: return(0);
223: }
225: PETSC_INTERN PetscErrorCode TaoCreate_BQNK(Tao tao)
226: {
227: TAO_BNK *bnk;
228: TAO_BQNK *bqnk;
232: TaoCreate_BNK(tao);
233: KSPSetOptionsPrefix(tao->ksp, "tao_bqnk_");
234: tao->ops->solve = TaoSolve_BQNK;
235: tao->ops->setfromoptions = TaoSetFromOptions_BQNK;
236: tao->ops->destroy = TaoDestroy_BQNK;
237: tao->ops->view = TaoView_BQNK;
238: tao->ops->setup = TaoSetUp_BQNK;
240: bnk = (TAO_BNK *)tao->data;
241: bnk->computehessian = TaoBQNKComputeHessian;
242: bnk->computestep = TaoBQNKComputeStep;
243: bnk->init_type = BNK_INIT_DIRECTION;
245: PetscNewLog(tao,&bqnk);
246: bnk->ctx = (void*)bqnk;
247: bqnk->is_spd = PETSC_TRUE;
249: MatCreate(PetscObjectComm((PetscObject)tao), &bqnk->B);
250: PetscObjectIncrementTabLevel((PetscObject)bqnk->B, (PetscObject)tao, 1);
251: MatSetOptionsPrefix(bqnk->B, "tao_bqnk_");
252: MatSetType(bqnk->B, MATLMVMSR1);
253: return(0);
254: }
256: /*@
257: TaoGetLMVMMatrix - Returns a pointer to the internal LMVM matrix. Valid
258: only for quasi-Newton family of methods.
260: Input Parameters:
261: . tao - Tao solver context
263: Output Parameters:
264: . B - LMVM matrix
266: Level: advanced
268: .seealso: TAOBQNLS, TAOBQNKLS, TAOBQNKTL, TAOBQNKTR, MATLMVM, TaoSetLMVMMatrix()
269: @*/
270: PetscErrorCode TaoGetLMVMMatrix(Tao tao, Mat *B)
271: {
272: TAO_BNK *bnk = (TAO_BNK*)tao->data;
273: TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx;
275: PetscBool flg = PETSC_FALSE;
278: PetscObjectTypeCompareAny((PetscObject)tao, &flg, TAOBQNLS, TAOBQNKLS, TAOBQNKTR, TAOBQNKTL, "");
279: if (!flg) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM Matrix only exists for quasi-Newton algorithms");
280: *B = bqnk->B;
281: return(0);
282: }
284: /*@
285: TaoSetLMVMMatrix - Sets an external LMVM matrix into the Tao solver. Valid
286: only for quasi-Newton family of methods.
288: QN family of methods create their own LMVM matrices and users who wish to
289: manipulate this matrix should use TaoGetLMVMMatrix() instead.
291: Input Parameters:
292: + tao - Tao solver context
293: - B - LMVM matrix
295: Level: advanced
297: .seealso: TAOBQNLS, TAOBQNKLS, TAOBQNKTL, TAOBQNKTR, MATLMVM, TaoGetLMVMMatrix()
298: @*/
299: PetscErrorCode TaoSetLMVMMatrix(Tao tao, Mat B)
300: {
301: TAO_BNK *bnk = (TAO_BNK*)tao->data;
302: TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx;
304: PetscBool flg = PETSC_FALSE;
307: PetscObjectTypeCompareAny((PetscObject)tao, &flg, TAOBQNLS, TAOBQNKLS, TAOBQNKTR, TAOBQNKTL, "");
308: if (!flg) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM Matrix only exists for quasi-Newton algorithms");
309: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &flg);
310: if (!flg) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "Given matrix is not an LMVM matrix");
311: if (bqnk->B) {
312: MatDestroy(&bqnk->B);
313: }
314: PetscObjectReference((PetscObject)B);
315: bqnk->B = B;
316: return(0);
317: }