1 #include <petsctaolinesearch.h> /*I "petsctaolinesearch.h" I*/ 2 #include <../src/tao/unconstrained/impls/lmvm/lmvm.h> 3 #include <../src/tao/bound/impls/blmvm/blmvm.h> 4 5 /*------------------------------------------------------------*/ 6 static PetscErrorCode TaoSolve_BLMVM(Tao tao) 7 { 8 TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data; 9 TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING; 10 PetscReal f, fold, gdx, gnorm, gnorm2; 11 PetscReal stepsize = 1.0,delta; 12 13 PetscFunctionBegin; 14 /* Project initial point onto bounds */ 15 PetscCall(TaoComputeVariableBounds(tao)); 16 PetscCall(VecMedian(tao->XL,tao->solution,tao->XU,tao->solution)); 17 PetscCall(TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU)); 18 19 /* Check convergence criteria */ 20 PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution,&f,blmP->unprojected_gradient)); 21 PetscCall(VecBoundGradientProjection(blmP->unprojected_gradient,tao->solution, tao->XL,tao->XU,tao->gradient)); 22 23 PetscCall(TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm)); 24 PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm),PetscObjectComm((PetscObject)tao),PETSC_ERR_USER, "User provided compute function generated Inf or NaN"); 25 26 tao->reason = TAO_CONTINUE_ITERATING; 27 PetscCall(TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its)); 28 PetscCall(TaoMonitor(tao,tao->niter,f,gnorm,0.0,stepsize)); 29 PetscCall((*tao->ops->convergencetest)(tao,tao->cnvP)); 30 if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 31 32 /* Set counter for gradient/reset steps */ 33 if (!blmP->recycle) { 34 blmP->grad = 0; 35 blmP->reset = 0; 36 PetscCall(MatLMVMReset(blmP->M, PETSC_FALSE)); 37 } 38 39 /* Have not converged; continue with Newton method */ 40 while (tao->reason == TAO_CONTINUE_ITERATING) { 41 /* Call general purpose update function */ 42 if (tao->ops->update) { 43 PetscCall((*tao->ops->update)(tao, tao->niter, tao->user_update)); 44 } 45 /* Compute direction */ 46 gnorm2 = gnorm*gnorm; 47 if (gnorm2 == 0.0) gnorm2 = PETSC_MACHINE_EPSILON; 48 if (f == 0.0) { 49 delta = 2.0 / gnorm2; 50 } else { 51 delta = 2.0 * PetscAbsScalar(f) / gnorm2; 52 } 53 PetscCall(MatLMVMSymBroydenSetDelta(blmP->M, delta)); 54 PetscCall(MatLMVMUpdate(blmP->M, tao->solution, tao->gradient)); 55 PetscCall(MatSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection)); 56 PetscCall(VecBoundGradientProjection(tao->stepdirection,tao->solution,tao->XL,tao->XU,tao->gradient)); 57 58 /* Check for success (descent direction) */ 59 PetscCall(VecDot(blmP->unprojected_gradient, tao->gradient, &gdx)); 60 if (gdx <= 0) { 61 /* Step is not descent or solve was not successful 62 Use steepest descent direction (scaled) */ 63 ++blmP->grad; 64 65 PetscCall(MatLMVMReset(blmP->M, PETSC_FALSE)); 66 PetscCall(MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient)); 67 PetscCall(MatSolve(blmP->M,blmP->unprojected_gradient, tao->stepdirection)); 68 } 69 PetscCall(VecScale(tao->stepdirection,-1.0)); 70 71 /* Perform the linesearch */ 72 fold = f; 73 PetscCall(VecCopy(tao->solution, blmP->Xold)); 74 PetscCall(VecCopy(blmP->unprojected_gradient, blmP->Gold)); 75 PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch,1.0)); 76 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, blmP->unprojected_gradient, tao->stepdirection, &stepsize, &ls_status)); 77 PetscCall(TaoAddLineSearchCounts(tao)); 78 79 if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) { 80 /* Linesearch failed 81 Reset factors and use scaled (projected) gradient step */ 82 ++blmP->reset; 83 84 f = fold; 85 PetscCall(VecCopy(blmP->Xold, tao->solution)); 86 PetscCall(VecCopy(blmP->Gold, blmP->unprojected_gradient)); 87 88 PetscCall(MatLMVMReset(blmP->M, PETSC_FALSE)); 89 PetscCall(MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient)); 90 PetscCall(MatSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection)); 91 PetscCall(VecScale(tao->stepdirection, -1.0)); 92 93 /* This may be incorrect; linesearch has values for stepmax and stepmin 94 that should be reset. */ 95 PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch,1.0)); 96 PetscCall(TaoLineSearchApply(tao->linesearch,tao->solution,&f, blmP->unprojected_gradient, tao->stepdirection, &stepsize, &ls_status)); 97 PetscCall(TaoAddLineSearchCounts(tao)); 98 99 if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) { 100 tao->reason = TAO_DIVERGED_LS_FAILURE; 101 break; 102 } 103 } 104 105 /* Check for converged */ 106 PetscCall(VecBoundGradientProjection(blmP->unprojected_gradient, tao->solution, tao->XL, tao->XU, tao->gradient)); 107 PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm)); 108 PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm),PetscObjectComm((PetscObject)tao),PETSC_ERR_USER, "User provided compute function generated Not-a-Number"); 109 tao->niter++; 110 PetscCall(TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its)); 111 PetscCall(TaoMonitor(tao,tao->niter,f,gnorm,0.0,stepsize)); 112 PetscCall((*tao->ops->convergencetest)(tao,tao->cnvP)); 113 } 114 PetscFunctionReturn(0); 115 } 116 117 static PetscErrorCode TaoSetup_BLMVM(Tao tao) 118 { 119 TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data; 120 121 PetscFunctionBegin; 122 /* Existence of tao->solution checked in TaoSetup() */ 123 PetscCall(VecDuplicate(tao->solution,&blmP->Xold)); 124 PetscCall(VecDuplicate(tao->solution,&blmP->Gold)); 125 PetscCall(VecDuplicate(tao->solution, &blmP->unprojected_gradient)); 126 if (!tao->stepdirection) { 127 PetscCall(VecDuplicate(tao->solution, &tao->stepdirection)); 128 } 129 if (!tao->gradient) { 130 PetscCall(VecDuplicate(tao->solution,&tao->gradient)); 131 } 132 /* Allocate matrix for the limited memory approximation */ 133 PetscCall(MatLMVMAllocate(blmP->M,tao->solution,blmP->unprojected_gradient)); 134 135 /* If the user has set a matrix to solve as the initial H0, set the options prefix here, and set up the KSP */ 136 if (blmP->H0) { 137 PetscCall(MatLMVMSetJ0(blmP->M, blmP->H0)); 138 } 139 PetscFunctionReturn(0); 140 } 141 142 /* ---------------------------------------------------------- */ 143 static PetscErrorCode TaoDestroy_BLMVM(Tao tao) 144 { 145 TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data; 146 147 PetscFunctionBegin; 148 if (tao->setupcalled) { 149 PetscCall(VecDestroy(&blmP->unprojected_gradient)); 150 PetscCall(VecDestroy(&blmP->Xold)); 151 PetscCall(VecDestroy(&blmP->Gold)); 152 } 153 PetscCall(MatDestroy(&blmP->M)); 154 if (blmP->H0) { 155 PetscObjectDereference((PetscObject)blmP->H0); 156 } 157 PetscCall(PetscFree(tao->data)); 158 PetscFunctionReturn(0); 159 } 160 161 /*------------------------------------------------------------*/ 162 static PetscErrorCode TaoSetFromOptions_BLMVM(PetscOptionItems* PetscOptionsObject,Tao tao) 163 { 164 TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data; 165 PetscBool is_spd; 166 167 PetscFunctionBegin; 168 PetscOptionsHeadBegin(PetscOptionsObject,"Limited-memory variable-metric method for bound constrained optimization"); 169 PetscCall(PetscOptionsBool("-tao_blmvm_recycle","enable recycling of the BFGS matrix between subsequent TaoSolve() calls","",blmP->recycle,&blmP->recycle,NULL)); 170 PetscOptionsHeadEnd(); 171 PetscCall(MatSetOptionsPrefix(blmP->M, ((PetscObject)tao)->prefix)); 172 PetscCall(MatAppendOptionsPrefix(blmP->M, "tao_blmvm_")); 173 PetscCall(MatSetFromOptions(blmP->M)); 174 PetscCall(MatGetOption(blmP->M, MAT_SPD, &is_spd)); 175 PetscCheck(is_spd,PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix must be symmetric positive-definite"); 176 PetscFunctionReturn(0); 177 } 178 179 /*------------------------------------------------------------*/ 180 static PetscErrorCode TaoView_BLMVM(Tao tao, PetscViewer viewer) 181 { 182 TAO_BLMVM *lmP = (TAO_BLMVM *)tao->data; 183 PetscBool isascii; 184 185 PetscFunctionBegin; 186 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 187 if (isascii) { 188 PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient steps: %" PetscInt_FMT "\n", lmP->grad)); 189 PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO)); 190 PetscCall(MatView(lmP->M, viewer)); 191 PetscCall(PetscViewerPopFormat(viewer)); 192 } 193 PetscFunctionReturn(0); 194 } 195 196 static PetscErrorCode TaoComputeDual_BLMVM(Tao tao, Vec DXL, Vec DXU) 197 { 198 TAO_BLMVM *blm = (TAO_BLMVM *) tao->data; 199 200 PetscFunctionBegin; 201 PetscValidHeaderSpecific(tao,TAO_CLASSID,1); 202 PetscValidHeaderSpecific(DXL,VEC_CLASSID,2); 203 PetscValidHeaderSpecific(DXU,VEC_CLASSID,3); 204 PetscCheck(tao->gradient && blm->unprojected_gradient,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Dual variables don't exist yet or no longer exist."); 205 206 PetscCall(VecCopy(tao->gradient,DXL)); 207 PetscCall(VecAXPY(DXL,-1.0,blm->unprojected_gradient)); 208 PetscCall(VecSet(DXU,0.0)); 209 PetscCall(VecPointwiseMax(DXL,DXL,DXU)); 210 211 PetscCall(VecCopy(blm->unprojected_gradient,DXU)); 212 PetscCall(VecAXPY(DXU,-1.0,tao->gradient)); 213 PetscCall(VecAXPY(DXU,1.0,DXL)); 214 PetscFunctionReturn(0); 215 } 216 217 /* ---------------------------------------------------------- */ 218 /*MC 219 TAOBLMVM - Bounded limited memory variable metric is a quasi-Newton method 220 for nonlinear minimization with bound constraints. It is an extension 221 of TAOLMVM 222 223 Options Database Keys: 224 . -tao_lmm_recycle - enable recycling of LMVM information between subsequent TaoSolve calls 225 226 Level: beginner 227 M*/ 228 PETSC_EXTERN PetscErrorCode TaoCreate_BLMVM(Tao tao) 229 { 230 TAO_BLMVM *blmP; 231 const char *morethuente_type = TAOLINESEARCHMT; 232 233 PetscFunctionBegin; 234 tao->ops->setup = TaoSetup_BLMVM; 235 tao->ops->solve = TaoSolve_BLMVM; 236 tao->ops->view = TaoView_BLMVM; 237 tao->ops->setfromoptions = TaoSetFromOptions_BLMVM; 238 tao->ops->destroy = TaoDestroy_BLMVM; 239 tao->ops->computedual = TaoComputeDual_BLMVM; 240 241 PetscCall(PetscNewLog(tao,&blmP)); 242 blmP->H0 = NULL; 243 blmP->recycle = PETSC_FALSE; 244 tao->data = (void*)blmP; 245 246 /* Override default settings (unless already changed) */ 247 if (!tao->max_it_changed) tao->max_it = 2000; 248 if (!tao->max_funcs_changed) tao->max_funcs = 4000; 249 250 PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch)); 251 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1)); 252 PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type)); 253 PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch,tao)); 254 255 PetscCall(KSPInitializePackage()); 256 PetscCall(MatCreate(((PetscObject)tao)->comm, &blmP->M)); 257 PetscCall(MatSetType(blmP->M, MATLMVMBFGS)); 258 PetscCall(PetscObjectIncrementTabLevel((PetscObject)blmP->M, (PetscObject)tao, 1)); 259 PetscFunctionReturn(0); 260 } 261 262 /*@ 263 TaoLMVMRecycle - Enable/disable recycling of the QN history between subsequent TaoSolve calls. 264 265 Input Parameters: 266 + tao - the Tao solver context 267 - flg - Boolean flag for recycling (PETSC_TRUE or PETSC_FALSE) 268 269 Level: intermediate 270 @*/ 271 PetscErrorCode TaoLMVMRecycle(Tao tao, PetscBool flg) 272 { 273 TAO_LMVM *lmP; 274 TAO_BLMVM *blmP; 275 PetscBool is_lmvm, is_blmvm; 276 277 PetscFunctionBegin; 278 PetscCall(PetscObjectTypeCompare((PetscObject)tao,TAOLMVM,&is_lmvm)); 279 PetscCall(PetscObjectTypeCompare((PetscObject)tao,TAOBLMVM,&is_blmvm)); 280 if (is_lmvm) { 281 lmP = (TAO_LMVM *)tao->data; 282 lmP->recycle = flg; 283 } else if (is_blmvm) { 284 blmP = (TAO_BLMVM *)tao->data; 285 blmP->recycle = flg; 286 } 287 PetscFunctionReturn(0); 288 } 289 290 /*@ 291 TaoLMVMSetH0 - Set the initial Hessian for the QN approximation 292 293 Input Parameters: 294 + tao - the Tao solver context 295 - H0 - Mat object for the initial Hessian 296 297 Level: advanced 298 299 .seealso: `TaoLMVMGetH0()`, `TaoLMVMGetH0KSP()` 300 @*/ 301 PetscErrorCode TaoLMVMSetH0(Tao tao, Mat H0) 302 { 303 TAO_LMVM *lmP; 304 TAO_BLMVM *blmP; 305 PetscBool is_lmvm, is_blmvm; 306 307 PetscFunctionBegin; 308 PetscCall(PetscObjectTypeCompare((PetscObject)tao,TAOLMVM,&is_lmvm)); 309 PetscCall(PetscObjectTypeCompare((PetscObject)tao,TAOBLMVM,&is_blmvm)); 310 if (is_lmvm) { 311 lmP = (TAO_LMVM *)tao->data; 312 PetscCall(PetscObjectReference((PetscObject)H0)); 313 lmP->H0 = H0; 314 } else if (is_blmvm) { 315 blmP = (TAO_BLMVM *)tao->data; 316 PetscCall(PetscObjectReference((PetscObject)H0)); 317 blmP->H0 = H0; 318 } 319 PetscFunctionReturn(0); 320 } 321 322 /*@ 323 TaoLMVMGetH0 - Get the matrix object for the QN initial Hessian 324 325 Input Parameters: 326 . tao - the Tao solver context 327 328 Output Parameters: 329 . H0 - Mat object for the initial Hessian 330 331 Level: advanced 332 333 .seealso: `TaoLMVMSetH0()`, `TaoLMVMGetH0KSP()` 334 @*/ 335 PetscErrorCode TaoLMVMGetH0(Tao tao, Mat *H0) 336 { 337 TAO_LMVM *lmP; 338 TAO_BLMVM *blmP; 339 PetscBool is_lmvm, is_blmvm; 340 Mat M; 341 342 PetscFunctionBegin; 343 PetscCall(PetscObjectTypeCompare((PetscObject)tao,TAOLMVM,&is_lmvm)); 344 PetscCall(PetscObjectTypeCompare((PetscObject)tao,TAOBLMVM,&is_blmvm)); 345 if (is_lmvm) { 346 lmP = (TAO_LMVM *)tao->data; 347 M = lmP->M; 348 } else if (is_blmvm) { 349 blmP = (TAO_BLMVM *)tao->data; 350 M = blmP->M; 351 } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONG, "This routine applies to TAO_LMVM and TAO_BLMVM."); 352 PetscCall(MatLMVMGetJ0(M, H0)); 353 PetscFunctionReturn(0); 354 } 355 356 /*@ 357 TaoLMVMGetH0KSP - Get the iterative solver for applying the inverse of the QN initial Hessian 358 359 Input Parameters: 360 . tao - the Tao solver context 361 362 Output Parameters: 363 . ksp - KSP solver context for the initial Hessian 364 365 Level: advanced 366 367 .seealso: `TaoLMVMGetH0()`, `TaoLMVMGetH0KSP()` 368 @*/ 369 PetscErrorCode TaoLMVMGetH0KSP(Tao tao, KSP *ksp) 370 { 371 TAO_LMVM *lmP; 372 TAO_BLMVM *blmP; 373 PetscBool is_lmvm, is_blmvm; 374 Mat M; 375 376 PetscFunctionBegin; 377 PetscCall(PetscObjectTypeCompare((PetscObject)tao,TAOLMVM,&is_lmvm)); 378 PetscCall(PetscObjectTypeCompare((PetscObject)tao,TAOBLMVM,&is_blmvm)); 379 if (is_lmvm) { 380 lmP = (TAO_LMVM *)tao->data; 381 M = lmP->M; 382 } else if (is_blmvm) { 383 blmP = (TAO_BLMVM *)tao->data; 384 M = blmP->M; 385 } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONG, "This routine applies to TAO_LMVM and TAO_BLMVM."); 386 PetscCall(MatLMVMGetJ0KSP(M, ksp)); 387 PetscFunctionReturn(0); 388 } 389