1 #include <../src/tao/bound/impls/bqnk/bqnk.h> 2 #include <petscksp.h> 3 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"}; 7 8 static PetscErrorCode TaoBQNKComputeHessian(Tao tao) 9 { 10 TAO_BNK *bnk = (TAO_BNK *)tao->data; 11 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 12 PetscErrorCode ierr; 13 PetscReal gnorm2, delta; 14 15 PetscFunctionBegin; 16 /* Alias the LMVM matrix into the TAO hessian */ 17 if (tao->hessian) { 18 ierr = MatDestroy(&tao->hessian);CHKERRQ(ierr); 19 } 20 if (tao->hessian_pre) { 21 ierr = MatDestroy(&tao->hessian_pre);CHKERRQ(ierr); 22 } 23 ierr = PetscObjectReference((PetscObject)bqnk->B);CHKERRQ(ierr); 24 tao->hessian = bqnk->B; 25 ierr = PetscObjectReference((PetscObject)bqnk->B);CHKERRQ(ierr); 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 ierr = MatSymBrdnSetDelta(bqnk->B, delta);CHKERRQ(ierr); 37 } 38 ierr = MatLMVMUpdate(tao->hessian, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 39 ierr = MatLMVMResetShift(tao->hessian);CHKERRQ(ierr); 40 /* Prepare the reduced sub-matrices for the inactive set */ 41 ierr = MatDestroy(&bnk->H_inactive);CHKERRQ(ierr); 42 if (bnk->active_idx) { 43 ierr = MatCreateSubMatrixVirtual(tao->hessian, bnk->inactive_idx, bnk->inactive_idx, &bnk->H_inactive);CHKERRQ(ierr); 44 ierr = PCLMVMSetIS(bqnk->pc, bnk->inactive_idx);CHKERRQ(ierr); 45 } else { 46 ierr = PetscObjectReference((PetscObject)tao->hessian);CHKERRQ(ierr); 47 bnk->H_inactive = tao->hessian; 48 ierr = PCLMVMClearIS(bqnk->pc);CHKERRQ(ierr); 49 } 50 ierr = MatDestroy(&bnk->Hpre_inactive);CHKERRQ(ierr); 51 ierr = PetscObjectReference((PetscObject)bnk->H_inactive);CHKERRQ(ierr); 52 bnk->Hpre_inactive = bnk->H_inactive; 53 PetscFunctionReturn(0); 54 } 55 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; 60 PetscErrorCode ierr; 61 62 PetscFunctionBegin; 63 ierr = TaoBNKComputeStep(tao, shift, ksp_reason, step_type);CHKERRQ(ierr); 64 if (*ksp_reason < 0) { 65 /* Krylov solver failed to converge so reset the LMVM matrix */ 66 ierr = MatLMVMReset(bqnk->B, PETSC_FALSE);CHKERRQ(ierr); 67 ierr = MatLMVMUpdate(bqnk->B, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 68 } 69 PetscFunctionReturn(0); 70 } 71 72 PetscErrorCode TaoSetUp_BQNK(Tao tao) 73 { 74 TAO_BNK *bnk = (TAO_BNK *)tao->data; 75 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 76 PetscErrorCode ierr; 77 PetscInt n, N; 78 PetscBool is_lmvm, is_sym, is_spd; 79 80 PetscFunctionBegin; 81 ierr = TaoSetUp_BNK(tao);CHKERRQ(ierr); 82 ierr = VecGetLocalSize(tao->solution,&n);CHKERRQ(ierr); 83 ierr = VecGetSize(tao->solution,&N);CHKERRQ(ierr); 84 ierr = MatSetSizes(bqnk->B, n, n, N, N);CHKERRQ(ierr); 85 ierr = MatLMVMAllocate(bqnk->B,tao->solution,bnk->unprojected_gradient);CHKERRQ(ierr); 86 ierr = PetscObjectBaseTypeCompare((PetscObject)bqnk->B, MATLMVM, &is_lmvm);CHKERRQ(ierr); 87 if (!is_lmvm) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "Matrix must be an LMVM-type"); 88 ierr = MatGetOption(bqnk->B, MAT_SYMMETRIC, &is_sym);CHKERRQ(ierr); 89 if (!is_sym) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix must be symmetric"); 90 ierr = MatGetOption(bqnk->B, MAT_SPD, &is_spd);CHKERRQ(ierr); 91 ierr = KSPGetPC(tao->ksp, &bqnk->pc);CHKERRQ(ierr); 92 ierr = PCSetType(bqnk->pc, PCLMVM);CHKERRQ(ierr); 93 ierr = PCLMVMSetMatLMVM(bqnk->pc, bqnk->B);CHKERRQ(ierr); 94 PetscFunctionReturn(0); 95 } 96 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; 101 PetscErrorCode ierr; 102 KSPType ksp_type; 103 104 PetscFunctionBegin; 105 ierr = PetscOptionsHead(PetscOptionsObject,"Quasi-Newton-Krylov method for bound constrained optimization");CHKERRQ(ierr); 106 ierr = PetscOptionsEList("-tao_bqnk_init_type", "radius initialization type", "", BQNK_INIT, BQNK_INIT_TYPES, BQNK_INIT[bnk->init_type], &bnk->init_type, 0);CHKERRQ(ierr); 107 ierr = PetscOptionsEList("-tao_bqnk_update_type", "radius update type", "", BNK_UPDATE, BNK_UPDATE_TYPES, BNK_UPDATE[bnk->update_type], &bnk->update_type, 0);CHKERRQ(ierr); 108 ierr = PetscOptionsEList("-tao_bqnk_as_type", "active set estimation method", "", BNK_AS, BNK_AS_TYPES, BNK_AS[bnk->as_type], &bnk->as_type, 0);CHKERRQ(ierr); 109 ierr = PetscOptionsReal("-tao_bqnk_sval", "(developer) Hessian perturbation starting value", "", bnk->sval, &bnk->sval,NULL);CHKERRQ(ierr); 110 ierr = PetscOptionsReal("-tao_bqnk_imin", "(developer) minimum initial Hessian perturbation", "", bnk->imin, &bnk->imin,NULL);CHKERRQ(ierr); 111 ierr = PetscOptionsReal("-tao_bqnk_imax", "(developer) maximum initial Hessian perturbation", "", bnk->imax, &bnk->imax,NULL);CHKERRQ(ierr); 112 ierr = PetscOptionsReal("-tao_bqnk_imfac", "(developer) initial merit factor for Hessian perturbation", "", bnk->imfac, &bnk->imfac,NULL);CHKERRQ(ierr); 113 ierr = PetscOptionsReal("-tao_bqnk_pmin", "(developer) minimum Hessian perturbation", "", bnk->pmin, &bnk->pmin,NULL);CHKERRQ(ierr); 114 ierr = PetscOptionsReal("-tao_bqnk_pmax", "(developer) maximum Hessian perturbation", "", bnk->pmax, &bnk->pmax,NULL);CHKERRQ(ierr); 115 ierr = PetscOptionsReal("-tao_bqnk_pgfac", "(developer) Hessian perturbation growth factor", "", bnk->pgfac, &bnk->pgfac,NULL);CHKERRQ(ierr); 116 ierr = PetscOptionsReal("-tao_bqnk_psfac", "(developer) Hessian perturbation shrink factor", "", bnk->psfac, &bnk->psfac,NULL);CHKERRQ(ierr); 117 ierr = PetscOptionsReal("-tao_bqnk_pmgfac", "(developer) merit growth factor for Hessian perturbation", "", bnk->pmgfac, &bnk->pmgfac,NULL);CHKERRQ(ierr); 118 ierr = PetscOptionsReal("-tao_bqnk_pmsfac", "(developer) merit shrink factor for Hessian perturbation", "", bnk->pmsfac, &bnk->pmsfac,NULL);CHKERRQ(ierr); 119 ierr = PetscOptionsReal("-tao_bqnk_eta1", "(developer) threshold for rejecting step (-tao_bqnk_update_type reduction)", "", bnk->eta1, &bnk->eta1,NULL);CHKERRQ(ierr); 120 ierr = PetscOptionsReal("-tao_bqnk_eta2", "(developer) threshold for accepting marginal step (-tao_bqnk_update_type reduction)", "", bnk->eta2, &bnk->eta2,NULL);CHKERRQ(ierr); 121 ierr = PetscOptionsReal("-tao_bqnk_eta3", "(developer) threshold for accepting reasonable step (-tao_bqnk_update_type reduction)", "", bnk->eta3, &bnk->eta3,NULL);CHKERRQ(ierr); 122 ierr = PetscOptionsReal("-tao_bqnk_eta4", "(developer) threshold for accepting good step (-tao_bqnk_update_type reduction)", "", bnk->eta4, &bnk->eta4,NULL);CHKERRQ(ierr); 123 ierr = PetscOptionsReal("-tao_bqnk_alpha1", "(developer) radius reduction factor for rejected step (-tao_bqnk_update_type reduction)", "", bnk->alpha1, &bnk->alpha1,NULL);CHKERRQ(ierr); 124 ierr = PetscOptionsReal("-tao_bqnk_alpha2", "(developer) radius reduction factor for marginally accepted bad step (-tao_bqnk_update_type reduction)", "", bnk->alpha2, &bnk->alpha2,NULL);CHKERRQ(ierr); 125 ierr = PetscOptionsReal("-tao_bqnk_alpha3", "(developer) radius increase factor for reasonable accepted step (-tao_bqnk_update_type reduction)", "", bnk->alpha3, &bnk->alpha3,NULL);CHKERRQ(ierr); 126 ierr = PetscOptionsReal("-tao_bqnk_alpha4", "(developer) radius increase factor for good accepted step (-tao_bqnk_update_type reduction)", "", bnk->alpha4, &bnk->alpha4,NULL);CHKERRQ(ierr); 127 ierr = PetscOptionsReal("-tao_bqnk_alpha5", "(developer) radius increase factor for very good accepted step (-tao_bqnk_update_type reduction)", "", bnk->alpha5, &bnk->alpha5,NULL);CHKERRQ(ierr); 128 ierr = PetscOptionsReal("-tao_bqnk_nu1", "(developer) threshold for small line-search step length (-tao_bqnk_update_type step)", "", bnk->nu1, &bnk->nu1,NULL);CHKERRQ(ierr); 129 ierr = PetscOptionsReal("-tao_bqnk_nu2", "(developer) threshold for reasonable line-search step length (-tao_bqnk_update_type step)", "", bnk->nu2, &bnk->nu2,NULL);CHKERRQ(ierr); 130 ierr = PetscOptionsReal("-tao_bqnk_nu3", "(developer) threshold for large line-search step length (-tao_bqnk_update_type step)", "", bnk->nu3, &bnk->nu3,NULL);CHKERRQ(ierr); 131 ierr = PetscOptionsReal("-tao_bqnk_nu4", "(developer) threshold for very large line-search step length (-tao_bqnk_update_type step)", "", bnk->nu4, &bnk->nu4,NULL);CHKERRQ(ierr); 132 ierr = 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);CHKERRQ(ierr); 133 ierr = PetscOptionsReal("-tao_bqnk_omega2", "(developer) radius reduction factor for small line-search step length (-tao_bqnk_update_type step)", "", bnk->omega2, &bnk->omega2,NULL);CHKERRQ(ierr); 134 ierr = PetscOptionsReal("-tao_bqnk_omega3", "(developer) radius factor for decent line-search step length (-tao_bqnk_update_type step)", "", bnk->omega3, &bnk->omega3,NULL);CHKERRQ(ierr); 135 ierr = PetscOptionsReal("-tao_bqnk_omega4", "(developer) radius increase factor for large line-search step length (-tao_bqnk_update_type step)", "", bnk->omega4, &bnk->omega4,NULL);CHKERRQ(ierr); 136 ierr = 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);CHKERRQ(ierr); 137 ierr = PetscOptionsReal("-tao_bqnk_mu1", "(developer) threshold for accepting very good step (-tao_bqnk_update_type interpolation)", "", bnk->mu1, &bnk->mu1,NULL);CHKERRQ(ierr); 138 ierr = PetscOptionsReal("-tao_bqnk_mu2", "(developer) threshold for accepting good step (-tao_bqnk_update_type interpolation)", "", bnk->mu2, &bnk->mu2,NULL);CHKERRQ(ierr); 139 ierr = PetscOptionsReal("-tao_bqnk_gamma1", "(developer) radius reduction factor for rejected very bad step (-tao_bqnk_update_type interpolation)", "", bnk->gamma1, &bnk->gamma1,NULL);CHKERRQ(ierr); 140 ierr = PetscOptionsReal("-tao_bqnk_gamma2", "(developer) radius reduction factor for rejected bad step (-tao_bqnk_update_type interpolation)", "", bnk->gamma2, &bnk->gamma2,NULL);CHKERRQ(ierr); 141 ierr = PetscOptionsReal("-tao_bqnk_gamma3", "(developer) radius increase factor for accepted good step (-tao_bqnk_update_type interpolation)", "", bnk->gamma3, &bnk->gamma3,NULL);CHKERRQ(ierr); 142 ierr = PetscOptionsReal("-tao_bqnk_gamma4", "(developer) radius increase factor for accepted very good step (-tao_bqnk_update_type interpolation)", "", bnk->gamma4, &bnk->gamma4,NULL);CHKERRQ(ierr); 143 ierr = PetscOptionsReal("-tao_bqnk_theta", "(developer) trust region interpolation factor (-tao_bqnk_update_type interpolation)", "", bnk->theta, &bnk->theta,NULL);CHKERRQ(ierr); 144 ierr = PetscOptionsReal("-tao_bqnk_min_radius", "(developer) lower bound on initial radius", "", bnk->min_radius, &bnk->min_radius,NULL);CHKERRQ(ierr); 145 ierr = PetscOptionsReal("-tao_bqnk_max_radius", "(developer) upper bound on radius", "", bnk->max_radius, &bnk->max_radius,NULL);CHKERRQ(ierr); 146 ierr = PetscOptionsReal("-tao_bqnk_epsilon", "(developer) tolerance used when computing actual and predicted reduction", "", bnk->epsilon, &bnk->epsilon,NULL);CHKERRQ(ierr); 147 ierr = PetscOptionsReal("-tao_bqnk_as_tol", "(developer) initial tolerance used when estimating actively bounded variables", "", bnk->as_tol, &bnk->as_tol,NULL);CHKERRQ(ierr); 148 ierr = PetscOptionsReal("-tao_bqnk_as_step", "(developer) step length used when estimating actively bounded variables", "", bnk->as_step, &bnk->as_step,NULL);CHKERRQ(ierr); 149 ierr = 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);CHKERRQ(ierr); 150 ierr = PetscOptionsTail();CHKERRQ(ierr); 151 ierr = TaoSetFromOptions(bnk->bncg);CHKERRQ(ierr); 152 ierr = TaoLineSearchSetFromOptions(tao->linesearch);CHKERRQ(ierr); 153 ierr = KSPSetFromOptions(tao->ksp);CHKERRQ(ierr); 154 ierr = KSPGetType(tao->ksp,&ksp_type);CHKERRQ(ierr); 155 ierr = PetscStrcmp(ksp_type,KSPCGNASH,&bnk->is_nash);CHKERRQ(ierr); 156 ierr = PetscStrcmp(ksp_type,KSPCGSTCG,&bnk->is_stcg);CHKERRQ(ierr); 157 ierr = PetscStrcmp(ksp_type,KSPCGGLTR,&bnk->is_gltr);CHKERRQ(ierr); 158 if (bnk->init_type == BNK_INIT_INTERPOLATION) bnk->init_type = BNK_INIT_DIRECTION; 159 ierr = MatSetFromOptions(bqnk->B);CHKERRQ(ierr); 160 ierr = MatGetOption(bqnk->B, MAT_SPD, &bqnk->is_spd);CHKERRQ(ierr); 161 PetscFunctionReturn(0); 162 } 163 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; 168 PetscErrorCode ierr; 169 PetscBool isascii; 170 171 PetscFunctionBegin; 172 ierr = TaoView_BNK(tao, viewer);CHKERRQ(ierr); 173 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 174 if (isascii) { 175 ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr); 176 ierr = MatView(bqnk->B, viewer);CHKERRQ(ierr); 177 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 178 } 179 PetscFunctionReturn(0); 180 } 181 182 static PetscErrorCode TaoDestroy_BQNK(Tao tao) 183 { 184 TAO_BNK *bnk = (TAO_BNK*)tao->data; 185 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 186 PetscErrorCode ierr; 187 188 PetscFunctionBegin; 189 ierr = MatDestroy(&bnk->Hpre_inactive);CHKERRQ(ierr); 190 ierr = MatDestroy(&bnk->H_inactive);CHKERRQ(ierr); 191 ierr = MatDestroy(&bqnk->B);CHKERRQ(ierr); 192 ierr = PetscFree(bnk->ctx);CHKERRQ(ierr); 193 ierr = TaoDestroy_BNK(tao);CHKERRQ(ierr); 194 PetscFunctionReturn(0); 195 } 196 197 PETSC_INTERN PetscErrorCode TaoCreate_BQNK(Tao tao) 198 { 199 TAO_BNK *bnk; 200 TAO_BQNK *bqnk; 201 PetscErrorCode ierr; 202 203 PetscFunctionBegin; 204 ierr = TaoCreate_BNK(tao);CHKERRQ(ierr); 205 ierr = KSPSetOptionsPrefix(tao->ksp, "tao_bqnk_");CHKERRQ(ierr); 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 ierr = PetscNewLog(tao,&bqnk);CHKERRQ(ierr); 217 bnk->ctx = (void*)bqnk; 218 bqnk->is_spd = PETSC_TRUE; 219 220 ierr = MatCreate(PetscObjectComm((PetscObject)tao), &bqnk->B);CHKERRQ(ierr); 221 ierr = PetscObjectIncrementTabLevel((PetscObject)bqnk->B, (PetscObject)tao, 1);CHKERRQ(ierr); 222 ierr = MatSetOptionsPrefix(bqnk->B, "tao_bqnk_");CHKERRQ(ierr); 223 ierr = MatSetType(bqnk->B, MATLMVMSR1);CHKERRQ(ierr); 224 PetscFunctionReturn(0); 225 } 226 227 PetscErrorCode TaoGetLMVMMatrix(Tao tao, Mat *B) 228 { 229 TAO_BNK *bnk = (TAO_BNK*)tao->data; 230 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 231 PetscErrorCode ierr; 232 PetscBool is_bqnls, is_bqnkls, is_bqnktr, is_bqnktl; 233 234 PetscFunctionBegin; 235 ierr = PetscObjectTypeCompare((PetscObject)tao, TAOBQNLS, &is_bqnls);CHKERRQ(ierr); 236 ierr = PetscObjectTypeCompare((PetscObject)tao, TAOBQNKLS, &is_bqnkls);CHKERRQ(ierr); 237 ierr = PetscObjectTypeCompare((PetscObject)tao, TAOBQNKTR, &is_bqnktr);CHKERRQ(ierr); 238 ierr = PetscObjectTypeCompare((PetscObject)tao, TAOBQNKTL, &is_bqnktl);CHKERRQ(ierr); 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 PetscFunctionReturn(0); 242 } 243