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