Actual source code: bqnk.c

petsc-3.11.3 2019-06-26
Report Typos and Errors
  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:     MatSymBrdnSetDelta(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 TaoSetUp_BQNK(Tao tao)
 73: {
 74:   TAO_BNK        *bnk = (TAO_BNK *)tao->data;
 75:   TAO_BQNK       *bqnk = (TAO_BQNK*)bnk->ctx;
 77:   PetscInt       n, N;
 78:   PetscBool      is_lmvm, is_sym, is_spd;

 81:   TaoSetUp_BNK(tao);
 82:   VecGetLocalSize(tao->solution,&n);
 83:   VecGetSize(tao->solution,&N);
 84:   MatSetSizes(bqnk->B, n, n, N, N);
 85:   MatLMVMAllocate(bqnk->B,tao->solution,bnk->unprojected_gradient);
 86:   PetscObjectBaseTypeCompare((PetscObject)bqnk->B, MATLMVM, &is_lmvm);
 87:   if (!is_lmvm) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "Matrix must be an LMVM-type");
 88:   MatGetOption(bqnk->B, MAT_SYMMETRIC, &is_sym);
 89:   if (!is_sym) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix must be symmetric");
 90:   MatGetOption(bqnk->B, MAT_SPD, &is_spd);
 91:   KSPGetPC(tao->ksp, &bqnk->pc);
 92:   PCSetType(bqnk->pc, PCLMVM);
 93:   PCLMVMSetMatLMVM(bqnk->pc, bqnk->B);
 94:   return(0);
 95: }

 97: static PetscErrorCode TaoSetFromOptions_BQNK(PetscOptionItems *PetscOptionsObject,Tao tao)
 98: {
 99:   TAO_BNK        *bnk = (TAO_BNK *)tao->data;
100:   TAO_BQNK       *bqnk = (TAO_BQNK*)bnk->ctx;
102:   KSPType        ksp_type;

105:   PetscOptionsHead(PetscOptionsObject,"Quasi-Newton-Krylov method for bound constrained optimization");
106:   PetscOptionsEList("-tao_bqnk_init_type", "radius initialization type", "", BQNK_INIT, BQNK_INIT_TYPES, BQNK_INIT[bnk->init_type], &bnk->init_type, 0);
107:   PetscOptionsEList("-tao_bqnk_update_type", "radius update type", "", BNK_UPDATE, BNK_UPDATE_TYPES, BNK_UPDATE[bnk->update_type], &bnk->update_type, 0);
108:   PetscOptionsEList("-tao_bqnk_as_type", "active set estimation method", "", BNK_AS, BNK_AS_TYPES, BNK_AS[bnk->as_type], &bnk->as_type, 0);
109:   PetscOptionsReal("-tao_bqnk_sval", "(developer) Hessian perturbation starting value", "", bnk->sval, &bnk->sval,NULL);
110:   PetscOptionsReal("-tao_bqnk_imin", "(developer) minimum initial Hessian perturbation", "", bnk->imin, &bnk->imin,NULL);
111:   PetscOptionsReal("-tao_bqnk_imax", "(developer) maximum initial Hessian perturbation", "", bnk->imax, &bnk->imax,NULL);
112:   PetscOptionsReal("-tao_bqnk_imfac", "(developer) initial merit factor for Hessian perturbation", "", bnk->imfac, &bnk->imfac,NULL);
113:   PetscOptionsReal("-tao_bqnk_pmin", "(developer) minimum Hessian perturbation", "", bnk->pmin, &bnk->pmin,NULL);
114:   PetscOptionsReal("-tao_bqnk_pmax", "(developer) maximum Hessian perturbation", "", bnk->pmax, &bnk->pmax,NULL);
115:   PetscOptionsReal("-tao_bqnk_pgfac", "(developer) Hessian perturbation growth factor", "", bnk->pgfac, &bnk->pgfac,NULL);
116:   PetscOptionsReal("-tao_bqnk_psfac", "(developer) Hessian perturbation shrink factor", "", bnk->psfac, &bnk->psfac,NULL);
117:   PetscOptionsReal("-tao_bqnk_pmgfac", "(developer) merit growth factor for Hessian perturbation", "", bnk->pmgfac, &bnk->pmgfac,NULL);
118:   PetscOptionsReal("-tao_bqnk_pmsfac", "(developer) merit shrink factor for Hessian perturbation", "", bnk->pmsfac, &bnk->pmsfac,NULL);
119:   PetscOptionsReal("-tao_bqnk_eta1", "(developer) threshold for rejecting step (-tao_bqnk_update_type reduction)", "", bnk->eta1, &bnk->eta1,NULL);
120:   PetscOptionsReal("-tao_bqnk_eta2", "(developer) threshold for accepting marginal step (-tao_bqnk_update_type reduction)", "", bnk->eta2, &bnk->eta2,NULL);
121:   PetscOptionsReal("-tao_bqnk_eta3", "(developer) threshold for accepting reasonable step (-tao_bqnk_update_type reduction)", "", bnk->eta3, &bnk->eta3,NULL);
122:   PetscOptionsReal("-tao_bqnk_eta4", "(developer) threshold for accepting good step (-tao_bqnk_update_type reduction)", "", bnk->eta4, &bnk->eta4,NULL);
123:   PetscOptionsReal("-tao_bqnk_alpha1", "(developer) radius reduction factor for rejected step (-tao_bqnk_update_type reduction)", "", bnk->alpha1, &bnk->alpha1,NULL);
124:   PetscOptionsReal("-tao_bqnk_alpha2", "(developer) radius reduction factor for marginally accepted bad step (-tao_bqnk_update_type reduction)", "", bnk->alpha2, &bnk->alpha2,NULL);
125:   PetscOptionsReal("-tao_bqnk_alpha3", "(developer) radius increase factor for reasonable accepted step (-tao_bqnk_update_type reduction)", "", bnk->alpha3, &bnk->alpha3,NULL);
126:   PetscOptionsReal("-tao_bqnk_alpha4", "(developer) radius increase factor for good accepted step (-tao_bqnk_update_type reduction)", "", bnk->alpha4, &bnk->alpha4,NULL);
127:   PetscOptionsReal("-tao_bqnk_alpha5", "(developer) radius increase factor for very good accepted step (-tao_bqnk_update_type reduction)", "", bnk->alpha5, &bnk->alpha5,NULL);
128:   PetscOptionsReal("-tao_bqnk_nu1", "(developer) threshold for small line-search step length (-tao_bqnk_update_type step)", "", bnk->nu1, &bnk->nu1,NULL);
129:   PetscOptionsReal("-tao_bqnk_nu2", "(developer) threshold for reasonable line-search step length (-tao_bqnk_update_type step)", "", bnk->nu2, &bnk->nu2,NULL);
130:   PetscOptionsReal("-tao_bqnk_nu3", "(developer) threshold for large line-search step length (-tao_bqnk_update_type step)", "", bnk->nu3, &bnk->nu3,NULL);
131:   PetscOptionsReal("-tao_bqnk_nu4", "(developer) threshold for very large line-search step length (-tao_bqnk_update_type step)", "", bnk->nu4, &bnk->nu4,NULL);
132:   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);
133:   PetscOptionsReal("-tao_bqnk_omega2", "(developer) radius reduction factor for small line-search step length (-tao_bqnk_update_type step)", "", bnk->omega2, &bnk->omega2,NULL);
134:   PetscOptionsReal("-tao_bqnk_omega3", "(developer) radius factor for decent line-search step length (-tao_bqnk_update_type step)", "", bnk->omega3, &bnk->omega3,NULL);
135:   PetscOptionsReal("-tao_bqnk_omega4", "(developer) radius increase factor for large line-search step length (-tao_bqnk_update_type step)", "", bnk->omega4, &bnk->omega4,NULL);
136:   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);
137:   PetscOptionsReal("-tao_bqnk_mu1", "(developer) threshold for accepting very good step (-tao_bqnk_update_type interpolation)", "", bnk->mu1, &bnk->mu1,NULL);
138:   PetscOptionsReal("-tao_bqnk_mu2", "(developer) threshold for accepting good step (-tao_bqnk_update_type interpolation)", "", bnk->mu2, &bnk->mu2,NULL);
139:   PetscOptionsReal("-tao_bqnk_gamma1", "(developer) radius reduction factor for rejected very bad step (-tao_bqnk_update_type interpolation)", "", bnk->gamma1, &bnk->gamma1,NULL);
140:   PetscOptionsReal("-tao_bqnk_gamma2", "(developer) radius reduction factor for rejected bad step (-tao_bqnk_update_type interpolation)", "", bnk->gamma2, &bnk->gamma2,NULL);
141:   PetscOptionsReal("-tao_bqnk_gamma3", "(developer) radius increase factor for accepted good step (-tao_bqnk_update_type interpolation)", "", bnk->gamma3, &bnk->gamma3,NULL);
142:   PetscOptionsReal("-tao_bqnk_gamma4", "(developer) radius increase factor for accepted very good step (-tao_bqnk_update_type interpolation)", "", bnk->gamma4, &bnk->gamma4,NULL);
143:   PetscOptionsReal("-tao_bqnk_theta", "(developer) trust region interpolation factor (-tao_bqnk_update_type interpolation)", "", bnk->theta, &bnk->theta,NULL);
144:   PetscOptionsReal("-tao_bqnk_min_radius", "(developer) lower bound on initial radius", "", bnk->min_radius, &bnk->min_radius,NULL);
145:   PetscOptionsReal("-tao_bqnk_max_radius", "(developer) upper bound on radius", "", bnk->max_radius, &bnk->max_radius,NULL);
146:   PetscOptionsReal("-tao_bqnk_epsilon", "(developer) tolerance used when computing actual and predicted reduction", "", bnk->epsilon, &bnk->epsilon,NULL);
147:   PetscOptionsReal("-tao_bqnk_as_tol", "(developer) initial tolerance used when estimating actively bounded variables", "", bnk->as_tol, &bnk->as_tol,NULL);
148:   PetscOptionsReal("-tao_bqnk_as_step", "(developer) step length used when estimating actively bounded variables", "", bnk->as_step, &bnk->as_step,NULL);
149:   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);
150:   PetscOptionsTail();
151:   TaoSetFromOptions(bnk->bncg);
152:   TaoLineSearchSetFromOptions(tao->linesearch);
153:   KSPSetFromOptions(tao->ksp);
154:   KSPGetType(tao->ksp,&ksp_type);
155:   PetscStrcmp(ksp_type,KSPCGNASH,&bnk->is_nash);
156:   PetscStrcmp(ksp_type,KSPCGSTCG,&bnk->is_stcg);
157:   PetscStrcmp(ksp_type,KSPCGGLTR,&bnk->is_gltr);
158:   if (bnk->init_type == BNK_INIT_INTERPOLATION) bnk->init_type = BNK_INIT_DIRECTION;
159:   MatSetFromOptions(bqnk->B);
160:   MatGetOption(bqnk->B, MAT_SPD, &bqnk->is_spd);
161:   return(0);
162: }

164: static PetscErrorCode TaoView_BQNK(Tao tao, PetscViewer viewer)
165: {
166:   TAO_BNK        *bnk = (TAO_BNK*)tao->data;
167:   TAO_BQNK       *bqnk = (TAO_BQNK*)bnk->ctx;
169:   PetscBool      isascii;

172:   TaoView_BNK(tao, viewer);
173:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
174:   if (isascii) {
175:     PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO);
176:     MatView(bqnk->B, viewer);
177:     PetscViewerPopFormat(viewer);
178:   }
179:   return(0);
180: }

182: static PetscErrorCode TaoDestroy_BQNK(Tao tao)
183: {
184:   TAO_BNK        *bnk = (TAO_BNK*)tao->data;
185:   TAO_BQNK       *bqnk = (TAO_BQNK*)bnk->ctx;
187: 
189:   MatDestroy(&bnk->Hpre_inactive);
190:   MatDestroy(&bnk->H_inactive);
191:   MatDestroy(&bqnk->B);
192:   PetscFree(bnk->ctx);
193:   TaoDestroy_BNK(tao);
194:   return(0);
195: }

197: PETSC_INTERN PetscErrorCode TaoCreate_BQNK(Tao tao)
198: {
199:   TAO_BNK        *bnk;
200:   TAO_BQNK       *bqnk;
202: 
204:   TaoCreate_BNK(tao);
205:   KSPSetOptionsPrefix(tao->ksp, "tao_bqnk_");
206:   tao->ops->setfromoptions = TaoSetFromOptions_BQNK;
207:   tao->ops->destroy = TaoDestroy_BQNK;
208:   tao->ops->view = TaoView_BQNK;
209:   tao->ops->setup = TaoSetUp_BQNK;
210: 
211:   bnk = (TAO_BNK *)tao->data;
212:   bnk->computehessian = TaoBQNKComputeHessian;
213:   bnk->computestep = TaoBQNKComputeStep;
214:   bnk->init_type = BNK_INIT_DIRECTION;
215: 
216:   PetscNewLog(tao,&bqnk);
217:   bnk->ctx = (void*)bqnk;
218:   bqnk->is_spd = PETSC_TRUE;
219: 
220:   MatCreate(PetscObjectComm((PetscObject)tao), &bqnk->B);
221:   PetscObjectIncrementTabLevel((PetscObject)bqnk->B, (PetscObject)tao, 1);
222:   MatSetOptionsPrefix(bqnk->B, "tao_bqnk_");
223:   MatSetType(bqnk->B, MATLMVMSR1);
224:   return(0);
225: }

227: PetscErrorCode TaoGetLMVMMatrix(Tao tao, Mat *B)
228: {
229:   TAO_BNK        *bnk = (TAO_BNK*)tao->data;
230:   TAO_BQNK       *bqnk = (TAO_BQNK*)bnk->ctx;
232:   PetscBool      is_bqnls, is_bqnkls, is_bqnktr, is_bqnktl;
233: 
235:   PetscObjectTypeCompare((PetscObject)tao, TAOBQNLS, &is_bqnls);
236:   PetscObjectTypeCompare((PetscObject)tao, TAOBQNKLS, &is_bqnkls);
237:   PetscObjectTypeCompare((PetscObject)tao, TAOBQNKTR, &is_bqnktr);
238:   PetscObjectTypeCompare((PetscObject)tao, TAOBQNKTL, &is_bqnktl);
239:   if (!is_bqnls && !is_bqnkls && !is_bqnktr && is_bqnktl) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM Matrix only exists for quasi-Newton algorithms");
240:   *B = bqnk->B;
241:   return(0);
242: }