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