1 #include <petsctaolinesearch.h> 2 #include <../src/tao/unconstrained/impls/lmvm/lmvm.h> 3 #include <../src/tao/bound/impls/blmvm/blmvm.h> 4 5 #define BLMVM_STEP_BFGS 0 6 #define BLMVM_STEP_GRAD 1 7 8 #define BLMVM_AS_NONE 0 9 #define BLMVM_AS_BERTSEKAS 1 10 #define BLMVM_AS_SIZE 2 11 12 static const char *BLMVM_AS_TYPE[64] = {"none", "bertsekas"}; 13 14 PETSC_INTERN PetscErrorCode TaoBLMVMEstimateActiveSet(Tao tao, PetscInt asType) 15 { 16 PetscErrorCode ierr; 17 TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data; 18 19 PetscFunctionBegin; 20 if (!tao->bounded) PetscFunctionReturn(0); 21 switch (asType) { 22 case BLMVM_AS_NONE: 23 ierr = ISDestroy(&blmP->inactive_idx);CHKERRQ(ierr); 24 ierr = VecWhichInactive(tao->XL, tao->solution, blmP->unprojected_gradient, tao->XU, PETSC_TRUE, &blmP->inactive_idx);CHKERRQ(ierr); 25 ierr = ISDestroy(&blmP->active_idx);CHKERRQ(ierr); 26 ierr = ISComplementVec(blmP->inactive_idx, tao->solution, &blmP->active_idx);CHKERRQ(ierr); 27 break; 28 29 case BLMVM_AS_BERTSEKAS: 30 /* Use gradient descent to estimate the active set */ 31 ierr = VecCopy(blmP->unprojected_gradient, blmP->W);CHKERRQ(ierr); 32 ierr = VecScale(blmP->W, -1.0);CHKERRQ(ierr); 33 ierr = TaoEstimateActiveBounds(tao->solution, tao->XL, tao->XU, blmP->unprojected_gradient, blmP->W, blmP->work, blmP->as_step, &blmP->as_tol, 34 &blmP->active_lower, &blmP->active_upper, &blmP->active_fixed, &blmP->active_idx, &blmP->inactive_idx);CHKERRQ(ierr); 35 break; 36 37 default: 38 break; 39 } 40 PetscFunctionReturn(0); 41 } 42 43 PETSC_INTERN PetscErrorCode TaoBLMVMBoundStep(Tao tao, PetscInt asType, Vec step) 44 { 45 PetscErrorCode ierr; 46 TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data; 47 48 PetscFunctionBegin; 49 if (!tao->bounded) PetscFunctionReturn(0); 50 switch (asType) { 51 case BLMVM_AS_NONE: 52 ierr = VecISSet(step, blmP->active_idx, 0.0);CHKERRQ(ierr); 53 break; 54 55 case BLMVM_AS_BERTSEKAS: 56 ierr = TaoBoundStep(tao->solution, tao->XL, tao->XU, blmP->active_lower, blmP->active_upper, blmP->active_fixed, 1.0, step);CHKERRQ(ierr); 57 break; 58 59 default: 60 break; 61 } 62 PetscFunctionReturn(0); 63 } 64 65 /*------------------------------------------------------------*/ 66 PETSC_INTERN PetscErrorCode TaoSolve_BLMVM(Tao tao) 67 { 68 PetscErrorCode ierr; 69 TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data; 70 TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING; 71 PetscReal f, fold, gdx, gnorm, resnorm; 72 PetscReal stepsize = 1.0; 73 PetscInt nDiff, nupdates, stepType = BLMVM_STEP_GRAD; 74 75 PetscFunctionBegin; 76 /* Project initial point onto bounds */ 77 ierr = TaoComputeVariableBounds(tao);CHKERRQ(ierr); 78 ierr = TaoBoundSolution(tao->solution, tao->XL,tao->XU, 0.0, &nDiff, tao->solution);CHKERRQ(ierr); 79 if (tao->bounded) { 80 ierr = TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);CHKERRQ(ierr); 81 } 82 83 /* Check convergence criteria */ 84 ierr = TaoComputeObjectiveAndGradient(tao, tao->solution,&f,blmP->unprojected_gradient);CHKERRQ(ierr); 85 ierr = VecBoundGradientProjection(blmP->unprojected_gradient,tao->solution, tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr); 86 87 ierr = TaoGradientNorm(tao, blmP->unprojected_gradient,NORM_2,&gnorm);CHKERRQ(ierr); 88 if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN"); 89 90 tao->reason = TAO_CONTINUE_ITERATING; 91 ierr = VecFischer(tao->solution, blmP->unprojected_gradient, tao->XL, tao->XU, blmP->W);CHKERRQ(ierr); 92 ierr = VecNorm(blmP->W, NORM_2, &resnorm);CHKERRQ(ierr); 93 ierr = TaoLogConvergenceHistory(tao,f,resnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 94 ierr = TaoMonitor(tao,tao->niter,f,resnorm,0.0,stepsize);CHKERRQ(ierr); 95 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 96 if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 97 98 /* Set counter for gradient/reset steps */ 99 if (!blmP->recycle) { 100 blmP->bfgs = 0; 101 blmP->grad = 0; 102 ierr = MatLMVMReset(blmP->M, PETSC_FALSE);CHKERRQ(ierr); 103 } 104 105 /* Have not converged; continue with Newton method */ 106 while (tao->reason == TAO_CONTINUE_ITERATING) { 107 /* Estimate active set at the current iterate */ 108 ierr = TaoBLMVMEstimateActiveSet(tao, blmP->as_type);CHKERRQ(ierr); 109 110 /* Compute direction */ 111 if (blmP->H0) { 112 ierr = MatLMVMSetJ0(blmP->M, blmP->H0);CHKERRQ(ierr); 113 stepType = BLMVM_STEP_BFGS; 114 } else if (!blmP->no_scale) { 115 ierr = MatLMVMSetJ0(blmP->M, blmP->Mscale);CHKERRQ(ierr); 116 } 117 ierr = MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient);CHKERRQ(ierr); 118 ierr = VecCopy(blmP->unprojected_gradient, blmP->red_grad);CHKERRQ(ierr); 119 ierr = VecISSet(blmP->red_grad, blmP->active_idx, 0.0);CHKERRQ(ierr); 120 ierr = MatSolve(blmP->M, blmP->red_grad, tao->stepdirection);CHKERRQ(ierr); 121 ierr = VecScale(tao->stepdirection,-1.0);CHKERRQ(ierr); 122 ierr = TaoBLMVMBoundStep(tao, blmP->as_type, tao->stepdirection);CHKERRQ(ierr); 123 ierr = MatLMVMGetUpdateCount(blmP->M, &nupdates);CHKERRQ(ierr); 124 if (nupdates > 0) stepType = BLMVM_STEP_BFGS; 125 126 /* Check for success (descent direction) */ 127 ierr = VecDot(tao->stepdirection, blmP->red_grad, &gdx);CHKERRQ(ierr); 128 if (gdx >= 0 || PetscIsInfOrNanReal(gdx)) { 129 /* Step is not descent or solve was not successful 130 Use steepest descent direction (scaled) */ 131 ierr = MatLMVMReset(blmP->M, PETSC_FALSE);CHKERRQ(ierr); 132 ierr = MatLMVMClearJ0(blmP->M);CHKERRQ(ierr); 133 ierr = MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient);CHKERRQ(ierr); 134 ierr = MatSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection);CHKERRQ(ierr); 135 ierr = VecScale(tao->stepdirection,-1.0);CHKERRQ(ierr); 136 ierr = TaoBLMVMBoundStep(tao, blmP->as_type, tao->stepdirection);CHKERRQ(ierr); 137 stepType = BLMVM_STEP_GRAD; 138 } 139 140 /* Perform the linesearch */ 141 fold = f; 142 ierr = VecCopy(tao->solution, blmP->Xold);CHKERRQ(ierr); 143 ierr = VecCopy(blmP->unprojected_gradient, blmP->Gold);CHKERRQ(ierr); 144 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, blmP->unprojected_gradient, tao->stepdirection, &stepsize, &ls_status);CHKERRQ(ierr); 145 ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr); 146 147 if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER && stepType != BLMVM_STEP_GRAD) { 148 /* Linesearch failed 149 Reset factors and use scaled (projected) gradient step */ 150 f = fold; 151 ierr = VecCopy(blmP->Xold, tao->solution);CHKERRQ(ierr); 152 ierr = VecCopy(blmP->Gold, blmP->unprojected_gradient);CHKERRQ(ierr); 153 154 ierr = MatLMVMReset(blmP->M, PETSC_FALSE);CHKERRQ(ierr); 155 ierr = MatLMVMClearJ0(blmP->M);CHKERRQ(ierr); 156 ierr = MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient);CHKERRQ(ierr); 157 ierr = MatSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection);CHKERRQ(ierr); 158 ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr); 159 ierr = TaoBLMVMBoundStep(tao, blmP->as_type, tao->stepdirection);CHKERRQ(ierr); 160 stepType = BLMVM_STEP_GRAD; 161 162 ierr = TaoLineSearchApply(tao->linesearch,tao->solution,&f, blmP->unprojected_gradient, tao->stepdirection, &stepsize, &ls_status);CHKERRQ(ierr); 163 ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr); 164 } 165 166 if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) { 167 /* Line search failed on a gradient step, so just mark reason for divergence */ 168 f = fold; 169 ierr = VecCopy(blmP->Xold, tao->solution);CHKERRQ(ierr); 170 ierr = VecCopy(blmP->Gold, blmP->unprojected_gradient);CHKERRQ(ierr); 171 tao->reason = TAO_DIVERGED_LS_FAILURE; 172 } else { 173 /* LS found valid step, so tally the step type and compute projected gradient */ 174 switch (stepType) { 175 case BLMVM_STEP_BFGS: 176 ++blmP->bfgs; 177 break; 178 case BLMVM_STEP_GRAD: 179 ++blmP->grad; 180 break; 181 default: 182 break; 183 } 184 ierr = VecBoundGradientProjection(blmP->unprojected_gradient, tao->solution, tao->XL, tao->XU, tao->gradient);CHKERRQ(ierr); 185 ierr = TaoGradientNorm(tao, blmP->unprojected_gradient,NORM_2,&gnorm);CHKERRQ(ierr); 186 if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Not-a-Number"); 187 } 188 189 /* Check for converged */ 190 tao->niter++; 191 ierr = VecFischer(tao->solution, blmP->unprojected_gradient, tao->XL, tao->XU, blmP->W);CHKERRQ(ierr); 192 ierr = VecNorm(blmP->W, NORM_2, &resnorm);CHKERRQ(ierr); 193 ierr = TaoLogConvergenceHistory(tao,f,resnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 194 ierr = TaoMonitor(tao,tao->niter,f,resnorm,0.0,stepsize);CHKERRQ(ierr); 195 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 196 } 197 PetscFunctionReturn(0); 198 } 199 200 PETSC_INTERN PetscErrorCode TaoSetup_BLMVM(Tao tao) 201 { 202 TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data; 203 PetscInt n,N; 204 PetscErrorCode ierr; 205 PetscBool is_spd, is_symbrdn; 206 207 PetscFunctionBegin; 208 /* Existence of tao->solution checked in TaoSetup() */ 209 ierr = VecDuplicate(tao->solution,&blmP->Xold);CHKERRQ(ierr); 210 ierr = VecDuplicate(tao->solution,&blmP->Gold);CHKERRQ(ierr); 211 ierr = VecDuplicate(tao->solution, &blmP->unprojected_gradient);CHKERRQ(ierr); 212 ierr = VecDuplicate(tao->solution, &blmP->W);CHKERRQ(ierr); 213 ierr = VecDuplicate(tao->solution, &blmP->work);CHKERRQ(ierr); 214 ierr = VecDuplicate(tao->solution, &blmP->red_grad);CHKERRQ(ierr); 215 216 if (!tao->stepdirection) { 217 ierr = VecDuplicate(tao->solution, &tao->stepdirection);CHKERRQ(ierr); 218 } 219 if (!tao->gradient) { 220 ierr = VecDuplicate(tao->solution,&tao->gradient);CHKERRQ(ierr); 221 } 222 /* Create matrix for the limited memory approximation */ 223 ierr = VecGetLocalSize(tao->solution,&n);CHKERRQ(ierr); 224 ierr = VecGetSize(tao->solution,&N);CHKERRQ(ierr); 225 ierr = MatSetSizes(blmP->M, n, n, N, N);CHKERRQ(ierr); 226 ierr = MatLMVMAllocate(blmP->M,tao->solution,blmP->unprojected_gradient);CHKERRQ(ierr); 227 ierr = MatGetOption(blmP->M, MAT_SPD, &is_spd);CHKERRQ(ierr); 228 if (!is_spd) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix is not symmetric positive-definite."); 229 ierr = PetscObjectTypeCompare((PetscObject)blmP->M, MATLMVMSYMBRDN, &is_symbrdn); 230 if (is_symbrdn) blmP->no_scale = PETSC_TRUE; /* makes no sense to scale L-SymBrdn with SymBrdn diagonal */ 231 232 /* If the user has set a matrix to solve as the initial H0, set the options prefix here, and set up the KSP */ 233 if (blmP->H0) { 234 ierr = MatLMVMSetJ0(blmP->M, blmP->H0);CHKERRQ(ierr); 235 } else if (!blmP->no_scale) { 236 if (!blmP->Mscale) { 237 ierr = MatCreateLMVMDiagBrdn(PetscObjectComm((PetscObject)blmP->M), n, N, &blmP->Mscale);CHKERRQ(ierr); 238 ierr = MatSetOptionsPrefix(blmP->Mscale, "tao_blmvm_scale_");CHKERRQ(ierr); 239 ierr = MatLMVMAllocate(blmP->Mscale, tao->solution, blmP->unprojected_gradient);CHKERRQ(ierr); 240 } 241 ierr = MatLMVMSetJ0(blmP->M, blmP->Mscale);CHKERRQ(ierr); 242 } 243 PetscFunctionReturn(0); 244 } 245 246 /* ---------------------------------------------------------- */ 247 PETSC_INTERN PetscErrorCode TaoDestroy_BLMVM(Tao tao) 248 { 249 TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data; 250 PetscErrorCode ierr; 251 252 PetscFunctionBegin; 253 if (tao->setupcalled) { 254 ierr = VecDestroy(&blmP->unprojected_gradient);CHKERRQ(ierr); 255 ierr = VecDestroy(&blmP->Xold);CHKERRQ(ierr); 256 ierr = VecDestroy(&blmP->Gold);CHKERRQ(ierr); 257 ierr = VecDestroy(&blmP->W);CHKERRQ(ierr); 258 ierr = VecDestroy(&blmP->work);CHKERRQ(ierr); 259 ierr = VecDestroy(&blmP->red_grad);CHKERRQ(ierr); 260 } 261 ierr = ISDestroy(&blmP->active_lower);CHKERRQ(ierr); 262 ierr = ISDestroy(&blmP->active_upper);CHKERRQ(ierr); 263 ierr = ISDestroy(&blmP->active_fixed);CHKERRQ(ierr); 264 ierr = ISDestroy(&blmP->active_idx);CHKERRQ(ierr); 265 ierr = ISDestroy(&blmP->inactive_idx);CHKERRQ(ierr); 266 ierr = MatDestroy(&blmP->M);CHKERRQ(ierr); 267 if (blmP->H0) { 268 PetscObjectDereference((PetscObject)blmP->H0); 269 } 270 if (blmP->Mscale) { 271 ierr = MatDestroy(&blmP->Mscale);CHKERRQ(ierr); 272 } 273 ierr = PetscFree(tao->data);CHKERRQ(ierr); 274 PetscFunctionReturn(0); 275 } 276 277 /*------------------------------------------------------------*/ 278 PETSC_INTERN PetscErrorCode TaoSetFromOptions_BLMVM(PetscOptionItems* PetscOptionsObject,Tao tao) 279 { 280 TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data; 281 PetscErrorCode ierr; 282 PetscBool is_lmvm, is_spd; 283 284 PetscFunctionBegin; 285 ierr = PetscOptionsHead(PetscOptionsObject,"Limited-memory variable-metric method for bound constrained optimization");CHKERRQ(ierr); 286 ierr = PetscOptionsBool("-tao_blmvm_recycle","enable recycling of the BFGS matrix between subsequent TaoSolve() calls","",blmP->recycle,&blmP->recycle,NULL);CHKERRQ(ierr); 287 ierr = PetscOptionsBool("-tao_blmvm_no_scale","(developer) disable the diagonal Broyden scaling of the BFGS approximation","",blmP->no_scale,&blmP->no_scale,NULL);CHKERRQ(ierr); 288 ierr = PetscOptionsEList("-tao_blmvm_as_type","active set estimation method", "", BLMVM_AS_TYPE, BLMVM_AS_SIZE, BLMVM_AS_TYPE[blmP->as_type], &blmP->as_type,NULL);CHKERRQ(ierr); 289 ierr = PetscOptionsReal("-tao_blmvm_as_tol", "initial tolerance used when estimating actively bounded variables","",blmP->as_tol,&blmP->as_tol,NULL);CHKERRQ(ierr); 290 ierr = PetscOptionsReal("-tao_blmvm_as_step", "step length used when estimating actively bounded variables","",blmP->as_step,&blmP->as_step,NULL);CHKERRQ(ierr); 291 ierr = PetscOptionsTail();CHKERRQ(ierr); 292 ierr = TaoLineSearchSetFromOptions(tao->linesearch);CHKERRQ(ierr); 293 ierr = MatSetFromOptions(blmP->M);CHKERRQ(ierr); 294 ierr = PetscObjectBaseTypeCompare((PetscObject)blmP->M, MATLMVM, &is_lmvm);CHKERRQ(ierr); 295 if (!is_lmvm) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "Matrix must be an LMVM-type"); 296 ierr = MatGetOption(blmP->M, MAT_SPD, &is_spd);CHKERRQ(ierr); 297 if (!is_spd) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix must be symmetric positive-definite"); 298 PetscFunctionReturn(0); 299 } 300 301 /*------------------------------------------------------------*/ 302 PETSC_INTERN PetscErrorCode TaoView_BLMVM(Tao tao, PetscViewer viewer) 303 { 304 TAO_BLMVM *lmP = (TAO_BLMVM *)tao->data; 305 PetscBool isascii; 306 PetscErrorCode ierr; 307 PetscInt recycled_its; 308 309 PetscFunctionBegin; 310 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr); 311 if (isascii) { 312 ierr = PetscViewerASCIIPrintf(viewer, " Gradient steps: %D\n", lmP->grad);CHKERRQ(ierr); 313 if (lmP->recycle) { 314 ierr = PetscViewerASCIIPrintf(viewer, " Recycle: on\n");CHKERRQ(ierr); 315 recycled_its = lmP->bfgs + lmP->grad; 316 ierr = PetscViewerASCIIPrintf(viewer, " Total recycled iterations: %D\n", recycled_its);CHKERRQ(ierr); 317 } 318 ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr); 319 ierr = MatView(lmP->M, viewer);CHKERRQ(ierr); 320 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 321 } 322 PetscFunctionReturn(0); 323 } 324 325 PETSC_INTERN PetscErrorCode TaoComputeDual_BLMVM(Tao tao, Vec DXL, Vec DXU) 326 { 327 TAO_BLMVM *blm = (TAO_BLMVM *) tao->data; 328 PetscErrorCode ierr; 329 330 PetscFunctionBegin; 331 PetscValidHeaderSpecific(tao,TAO_CLASSID,1); 332 PetscValidHeaderSpecific(DXL,VEC_CLASSID,2); 333 PetscValidHeaderSpecific(DXU,VEC_CLASSID,3); 334 if (!tao->gradient || !blm->unprojected_gradient) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Dual variables don't exist yet or no longer exist.\n"); 335 336 ierr = VecCopy(tao->gradient,DXL);CHKERRQ(ierr); 337 ierr = VecAXPY(DXL,-1.0,blm->unprojected_gradient);CHKERRQ(ierr); 338 ierr = VecSet(DXU,0.0);CHKERRQ(ierr); 339 ierr = VecPointwiseMax(DXL,DXL,DXU);CHKERRQ(ierr); 340 341 ierr = VecCopy(blm->unprojected_gradient,DXU);CHKERRQ(ierr); 342 ierr = VecAXPY(DXU,-1.0,tao->gradient);CHKERRQ(ierr); 343 ierr = VecAXPY(DXU,1.0,DXL);CHKERRQ(ierr); 344 PetscFunctionReturn(0); 345 } 346 347 /* ---------------------------------------------------------- */ 348 /*MC 349 TAOBLMVM - Bounded limited memory variable metric is a quasi-Newton method 350 for nonlinear minimization with bound constraints. It is an extension 351 of TAOLMVM 352 353 Options Database Keys: 354 . -tao_blmvm_recycle - enable recycling LMVM updates between TaoSolve() calls 355 . -tao_blmvm_no_scale - (developer) disables diagonal Broyden scaling on the LMVM approximation 356 . -tao_blmvm_as_type - (developer) active set estimation method <none,bertsekas> 357 . -tao_blmvm_as_tol - (developer) initial distance tolerance used when estimating the active set 358 . -tao_blmvm_as_step - (developer) trial step length used when estimating the active set 359 360 Level: beginner 361 M*/ 362 PETSC_EXTERN PetscErrorCode TaoCreate_BLMVM(Tao tao) 363 { 364 TAO_BLMVM *blmP; 365 const char *morethuente_type = TAOLINESEARCHMT; 366 PetscErrorCode ierr; 367 368 PetscFunctionBegin; 369 tao->ops->setup = TaoSetup_BLMVM; 370 tao->ops->solve = TaoSolve_BLMVM; 371 tao->ops->view = TaoView_BLMVM; 372 tao->ops->setfromoptions = TaoSetFromOptions_BLMVM; 373 tao->ops->destroy = TaoDestroy_BLMVM; 374 tao->ops->computedual = TaoComputeDual_BLMVM; 375 376 ierr = PetscNewLog(tao,&blmP);CHKERRQ(ierr); 377 blmP->H0 = NULL; 378 blmP->as_step = 0.001; 379 blmP->as_tol = 0.001; 380 blmP->as_type = BLMVM_AS_BERTSEKAS; 381 blmP->recycle = PETSC_FALSE; 382 blmP->no_scale = PETSC_FALSE; 383 384 tao->data = (void*)blmP; 385 386 /* Override default settings (unless already changed) */ 387 if (!tao->max_it_changed) tao->max_it = 2000; 388 if (!tao->max_funcs_changed) tao->max_funcs = 4000; 389 390 ierr = TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);CHKERRQ(ierr); 391 ierr = PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);CHKERRQ(ierr); 392 ierr = TaoLineSearchSetType(tao->linesearch, morethuente_type);CHKERRQ(ierr); 393 ierr = TaoLineSearchUseTaoRoutines(tao->linesearch,tao);CHKERRQ(ierr); 394 ierr = TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);CHKERRQ(ierr); 395 396 ierr = KSPInitializePackage();CHKERRQ(ierr); 397 ierr = MatCreate(((PetscObject)tao)->comm, &blmP->M);CHKERRQ(ierr); 398 ierr = PetscObjectIncrementTabLevel((PetscObject)blmP->M, (PetscObject)tao, 1);CHKERRQ(ierr); 399 ierr = MatSetType(blmP->M, MATLMVMBFGS);CHKERRQ(ierr); 400 ierr = MatSetOptionsPrefix(blmP->M, "tao_blmvm_");CHKERRQ(ierr); 401 PetscFunctionReturn(0); 402 } 403 404 PETSC_EXTERN PetscErrorCode TaoBLMVMSetH0(Tao tao, Mat H0) 405 { 406 TAO_LMVM *lmP; 407 TAO_BLMVM *blmP; 408 TaoType type; 409 PetscBool is_lmvm, is_blmvm; 410 PetscErrorCode ierr; 411 412 ierr = TaoGetType(tao, &type);CHKERRQ(ierr); 413 ierr = PetscStrcmp(type, TAOLMVM, &is_lmvm);CHKERRQ(ierr); 414 ierr = PetscStrcmp(type, TAOBLMVM, &is_blmvm);CHKERRQ(ierr); 415 416 if (is_lmvm) { 417 lmP = (TAO_LMVM *)tao->data; 418 ierr = PetscObjectReference((PetscObject)H0);CHKERRQ(ierr); 419 lmP->H0 = H0; 420 } else if (is_blmvm) { 421 blmP = (TAO_BLMVM *)tao->data; 422 ierr = PetscObjectReference((PetscObject)H0);CHKERRQ(ierr); 423 blmP->H0 = H0; 424 } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONGSTATE, "This routine applies to TAO_LMVM and TAO_BLMVM."); 425 PetscFunctionReturn(0); 426 } 427 428 PETSC_EXTERN PetscErrorCode TaoBLMVMGetH0(Tao tao, Mat *H0) 429 { 430 TAO_LMVM *lmP; 431 TAO_BLMVM *blmP; 432 TaoType type; 433 PetscBool is_lmvm, is_blmvm; 434 Mat M; 435 436 PetscErrorCode ierr; 437 438 ierr = TaoGetType(tao, &type);CHKERRQ(ierr); 439 ierr = PetscStrcmp(type, TAOLMVM, &is_lmvm);CHKERRQ(ierr); 440 ierr = PetscStrcmp(type, TAOBLMVM, &is_blmvm);CHKERRQ(ierr); 441 442 if (is_lmvm) { 443 lmP = (TAO_LMVM *)tao->data; 444 M = lmP->M; 445 } else if (is_blmvm) { 446 blmP = (TAO_BLMVM *)tao->data; 447 M = blmP->M; 448 } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONGSTATE, "This routine applies to TAO_LMVM and TAO_BLMVM."); 449 ierr = MatLMVMGetJ0(M, H0);CHKERRQ(ierr); 450 PetscFunctionReturn(0); 451 } 452 453 PETSC_EXTERN PetscErrorCode TaoBLMVMGetH0KSP(Tao tao, KSP *ksp) 454 { 455 TAO_LMVM *lmP; 456 TAO_BLMVM *blmP; 457 TaoType type; 458 PetscBool is_lmvm, is_blmvm; 459 Mat M; 460 PetscErrorCode ierr; 461 462 ierr = TaoGetType(tao, &type);CHKERRQ(ierr); 463 ierr = PetscStrcmp(type, TAOLMVM, &is_lmvm);CHKERRQ(ierr); 464 ierr = PetscStrcmp(type, TAOBLMVM, &is_blmvm);CHKERRQ(ierr); 465 466 if (is_lmvm) { 467 lmP = (TAO_LMVM *)tao->data; 468 M = lmP->M; 469 } else if (is_blmvm) { 470 blmP = (TAO_BLMVM *)tao->data; 471 M = blmP->M; 472 } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONGSTATE, "This routine applies to TAO_LMVM and TAO_BLMVM."); 473 ierr = MatLMVMGetJ0KSP(M, ksp);CHKERRQ(ierr); 474 PetscFunctionReturn(0); 475 } 476