1 #include <petsctaolinesearch.h> 2 #include <../src/tao/matrix/lmvmmat.h> 3 #include <../src/tao/bound/impls/blmvm/blmvm.h> 4 5 /*------------------------------------------------------------*/ 6 #undef __FUNCT__ 7 #define __FUNCT__ "TaoSolve_BLMVM" 8 static PetscErrorCode TaoSolve_BLMVM(Tao tao) 9 { 10 PetscErrorCode ierr; 11 TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data; 12 TaoConvergedReason reason = TAO_CONTINUE_ITERATING; 13 TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING; 14 PetscReal f, fold, gdx, gnorm; 15 PetscReal stepsize = 1.0,delta; 16 PetscInt iter = 0; 17 18 PetscFunctionBegin; 19 /* Project initial point onto bounds */ 20 ierr = TaoComputeVariableBounds(tao);CHKERRQ(ierr); 21 ierr = VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);CHKERRQ(ierr); 22 ierr = TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);CHKERRQ(ierr); 23 24 /* Check convergence criteria */ 25 ierr = TaoComputeObjectiveAndGradient(tao, tao->solution,&f,blmP->unprojected_gradient);CHKERRQ(ierr); 26 ierr = VecBoundGradientProjection(blmP->unprojected_gradient,tao->solution, tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr); 27 28 ierr = VecNorm(tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr); 29 if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf pr NaN"); 30 31 ierr = TaoMonitor(tao, iter, f, gnorm, 0.0, stepsize, &reason);CHKERRQ(ierr); 32 if (reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 33 34 /* Set initial scaling for the function */ 35 if (f != 0.0) { 36 delta = 2.0*PetscAbsScalar(f) / (gnorm*gnorm); 37 } else { 38 delta = 2.0 / (gnorm*gnorm); 39 } 40 ierr = MatLMVMSetDelta(blmP->M,delta);CHKERRQ(ierr); 41 42 /* Set counter for gradient/reset steps */ 43 blmP->grad = 0; 44 blmP->reset = 0; 45 46 /* Have not converged; continue with Newton method */ 47 while (reason == TAO_CONTINUE_ITERATING) { 48 /* Compute direction */ 49 ierr = MatLMVMUpdate(blmP->M, tao->solution, tao->gradient);CHKERRQ(ierr); 50 ierr = MatLMVMSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection);CHKERRQ(ierr); 51 ierr = VecBoundGradientProjection(tao->stepdirection,tao->solution,tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr); 52 53 /* Check for success (descent direction) */ 54 ierr = VecDot(blmP->unprojected_gradient, tao->gradient, &gdx);CHKERRQ(ierr); 55 if (gdx <= 0) { 56 /* Step is not descent or solve was not successful 57 Use steepest descent direction (scaled) */ 58 ++blmP->grad; 59 60 if (f != 0.0) { 61 delta = 2.0*PetscAbsScalar(f) / (gnorm*gnorm); 62 } else { 63 delta = 2.0 / (gnorm*gnorm); 64 } 65 ierr = MatLMVMSetDelta(blmP->M,delta);CHKERRQ(ierr); 66 ierr = MatLMVMReset(blmP->M);CHKERRQ(ierr); 67 ierr = MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient);CHKERRQ(ierr); 68 ierr = MatLMVMSolve(blmP->M,blmP->unprojected_gradient, tao->stepdirection);CHKERRQ(ierr); 69 } 70 ierr = VecScale(tao->stepdirection,-1.0);CHKERRQ(ierr); 71 72 /* Perform the linesearch */ 73 fold = f; 74 ierr = VecCopy(tao->solution, blmP->Xold);CHKERRQ(ierr); 75 ierr = VecCopy(blmP->unprojected_gradient, blmP->Gold);CHKERRQ(ierr); 76 ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);CHKERRQ(ierr); 77 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, blmP->unprojected_gradient, tao->stepdirection, &stepsize, &ls_status);CHKERRQ(ierr); 78 ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr); 79 80 if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) { 81 /* Linesearch failed 82 Reset factors and use scaled (projected) gradient step */ 83 ++blmP->reset; 84 85 f = fold; 86 ierr = VecCopy(blmP->Xold, tao->solution);CHKERRQ(ierr); 87 ierr = VecCopy(blmP->Gold, blmP->unprojected_gradient);CHKERRQ(ierr); 88 89 if (f != 0.0) { 90 delta = 2.0* PetscAbsScalar(f) / (gnorm*gnorm); 91 } else { 92 delta = 2.0/ (gnorm*gnorm); 93 } 94 ierr = MatLMVMSetDelta(blmP->M,delta);CHKERRQ(ierr); 95 ierr = MatLMVMReset(blmP->M);CHKERRQ(ierr); 96 ierr = MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient);CHKERRQ(ierr); 97 ierr = MatLMVMSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection);CHKERRQ(ierr); 98 ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr); 99 100 /* This may be incorrect; linesearch has values fo stepmax and stepmin 101 that should be reset. */ 102 ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,1.0); 103 ierr = TaoLineSearchApply(tao->linesearch,tao->solution,&f, blmP->unprojected_gradient, tao->stepdirection, &stepsize, &ls_status);CHKERRQ(ierr); 104 ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr); 105 106 if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) { 107 tao->reason = TAO_DIVERGED_LS_FAILURE; 108 break; 109 } 110 } 111 112 /* Check for converged */ 113 ierr = VecBoundGradientProjection(blmP->unprojected_gradient, tao->solution, tao->XL, tao->XU, tao->gradient);CHKERRQ(ierr); 114 ierr = VecNorm(tao->gradient, NORM_2, &gnorm);CHKERRQ(ierr); 115 116 117 if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Not-a-Number"); 118 iter++; 119 ierr = TaoMonitor(tao, iter, f, gnorm, 0.0, stepsize, &reason);CHKERRQ(ierr); 120 } 121 PetscFunctionReturn(0); 122 } 123 124 #undef __FUNCT__ 125 #define __FUNCT__ "TaoSetup_BLMVM" 126 static PetscErrorCode TaoSetup_BLMVM(Tao tao) 127 { 128 TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data; 129 PetscInt n,N; 130 PetscErrorCode ierr; 131 132 PetscFunctionBegin; 133 /* Existence of tao->solution checked in TaoSetup() */ 134 ierr = VecDuplicate(tao->solution,&blmP->Xold);CHKERRQ(ierr); 135 ierr = VecDuplicate(tao->solution,&blmP->Gold);CHKERRQ(ierr); 136 ierr = VecDuplicate(tao->solution, &blmP->unprojected_gradient); 137 138 if (!tao->stepdirection) { 139 ierr = VecDuplicate(tao->solution, &tao->stepdirection);CHKERRQ(ierr); 140 } 141 if (!tao->gradient) { 142 ierr = VecDuplicate(tao->solution,&tao->gradient);CHKERRQ(ierr); 143 } 144 if (!tao->XL) { 145 ierr = VecDuplicate(tao->solution,&tao->XL);CHKERRQ(ierr); 146 ierr = VecSet(tao->XL,PETSC_NINFINITY);CHKERRQ(ierr); 147 } 148 if (!tao->XU) { 149 ierr = VecDuplicate(tao->solution,&tao->XU);CHKERRQ(ierr); 150 ierr = VecSet(tao->XU,PETSC_INFINITY);CHKERRQ(ierr); 151 } 152 /* Create matrix for the limited memory approximation */ 153 ierr = VecGetLocalSize(tao->solution,&n);CHKERRQ(ierr); 154 ierr = VecGetSize(tao->solution,&N);CHKERRQ(ierr); 155 ierr = MatCreateLMVM(((PetscObject)tao)->comm,n,N,&blmP->M);CHKERRQ(ierr); 156 ierr = MatLMVMAllocateVectors(blmP->M,tao->solution);CHKERRQ(ierr); 157 PetscFunctionReturn(0); 158 } 159 160 /* ---------------------------------------------------------- */ 161 #undef __FUNCT__ 162 #define __FUNCT__ "TaoDestroy_BLMVM" 163 static PetscErrorCode TaoDestroy_BLMVM(Tao tao) 164 { 165 TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data; 166 PetscErrorCode ierr; 167 168 PetscFunctionBegin; 169 if (tao->setupcalled) { 170 ierr = MatDestroy(&blmP->M);CHKERRQ(ierr); 171 ierr = VecDestroy(&blmP->unprojected_gradient);CHKERRQ(ierr); 172 ierr = VecDestroy(&blmP->Xold);CHKERRQ(ierr); 173 ierr = VecDestroy(&blmP->Gold);CHKERRQ(ierr); 174 } 175 ierr = PetscFree(tao->data);CHKERRQ(ierr); 176 PetscFunctionReturn(0); 177 } 178 179 /*------------------------------------------------------------*/ 180 #undef __FUNCT__ 181 #define __FUNCT__ "TaoSetFromOptions_BLMVM" 182 static PetscErrorCode TaoSetFromOptions_BLMVM(Tao tao) 183 { 184 PetscErrorCode ierr; 185 186 PetscFunctionBegin; 187 ierr = PetscOptionsHead("Limited-memory variable-metric method for bound constrained optimization");CHKERRQ(ierr); 188 ierr = TaoLineSearchSetFromOptions(tao->linesearch);CHKERRQ(ierr); 189 ierr = PetscOptionsTail();CHKERRQ(ierr); 190 PetscFunctionReturn(0); 191 } 192 193 194 /*------------------------------------------------------------*/ 195 #undef __FUNCT__ 196 #define __FUNCT__ "TaoView_BLMVM" 197 static int TaoView_BLMVM(Tao tao, PetscViewer viewer) 198 { 199 TAO_BLMVM *lmP = (TAO_BLMVM *)tao->data; 200 PetscBool isascii; 201 PetscErrorCode ierr; 202 203 PetscFunctionBegin; 204 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr); 205 if (isascii) { 206 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 207 ierr = PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", lmP->grad);CHKERRQ(ierr); 208 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 209 } 210 PetscFunctionReturn(0); 211 } 212 213 #undef __FUNCT__ 214 #define __FUNCT__ "TaoComputeDual_BLMVM" 215 static PetscErrorCode TaoComputeDual_BLMVM(Tao tao, Vec DXL, Vec DXU) 216 { 217 TAO_BLMVM *blm = (TAO_BLMVM *) tao->data; 218 PetscErrorCode ierr; 219 220 PetscFunctionBegin; 221 PetscValidHeaderSpecific(tao,TAO_CLASSID,1); 222 PetscValidHeaderSpecific(DXL,VEC_CLASSID,2); 223 PetscValidHeaderSpecific(DXU,VEC_CLASSID,3); 224 if (!tao->gradient || !blm->unprojected_gradient) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Dual variables don't exist yet or no longer exist.\n"); 225 226 ierr = VecCopy(tao->gradient,DXL);CHKERRQ(ierr); 227 ierr = VecAXPY(DXL,-1.0,blm->unprojected_gradient);CHKERRQ(ierr); 228 ierr = VecSet(DXU,0.0);CHKERRQ(ierr); 229 ierr = VecPointwiseMax(DXL,DXL,DXU);CHKERRQ(ierr); 230 231 ierr = VecCopy(blm->unprojected_gradient,DXU);CHKERRQ(ierr); 232 ierr = VecAXPY(DXU,-1.0,tao->gradient);CHKERRQ(ierr); 233 ierr = VecAXPY(DXU,1.0,DXL);CHKERRQ(ierr); 234 PetscFunctionReturn(0); 235 } 236 237 /* ---------------------------------------------------------- */ 238 /*MC 239 TAOBLMVM - Bounded limited memory variable metric is a quasi-Newton method 240 for nonlinear minimization with bound constraints. It is an extension 241 of TAOLMVM 242 243 Options Database Keys: 244 + -tao_lmm_vectors - number of vectors to use for approximation 245 . -tao_lmm_scale_type - "none","scalar","broyden" 246 . -tao_lmm_limit_type - "none","average","relative","absolute" 247 . -tao_lmm_rescale_type - "none","scalar","gl" 248 . -tao_lmm_limit_mu - mu limiting factor 249 . -tao_lmm_limit_nu - nu limiting factor 250 . -tao_lmm_delta_min - minimum delta value 251 . -tao_lmm_delta_max - maximum delta value 252 . -tao_lmm_broyden_phi - phi factor for Broyden scaling 253 . -tao_lmm_scalar_alpha - alpha factor for scalar scaling 254 . -tao_lmm_rescale_alpha - alpha factor for rescaling diagonal 255 . -tao_lmm_rescale_beta - beta factor for rescaling diagonal 256 . -tao_lmm_scalar_history - amount of history for scalar scaling 257 . -tao_lmm_rescale_history - amount of history for rescaling diagonal 258 - -tao_lmm_eps - rejection tolerance 259 260 Level: beginner 261 M*/ 262 #undef __FUNCT__ 263 #define __FUNCT__ "TaoCreate_BLMVM" 264 PETSC_EXTERN PetscErrorCode TaoCreate_BLMVM(Tao tao) 265 { 266 TAO_BLMVM *blmP; 267 const char *morethuente_type = TAOLINESEARCHMT; 268 PetscErrorCode ierr; 269 270 PetscFunctionBegin; 271 tao->ops->setup = TaoSetup_BLMVM; 272 tao->ops->solve = TaoSolve_BLMVM; 273 tao->ops->view = TaoView_BLMVM; 274 tao->ops->setfromoptions = TaoSetFromOptions_BLMVM; 275 tao->ops->destroy = TaoDestroy_BLMVM; 276 tao->ops->computedual = TaoComputeDual_BLMVM; 277 278 ierr = PetscNewLog(tao,&blmP);CHKERRQ(ierr); 279 tao->data = (void*)blmP; 280 tao->max_it = 2000; 281 tao->max_funcs = 4000; 282 tao->fatol = 1e-4; 283 tao->frtol = 1e-4; 284 285 ierr = TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);CHKERRQ(ierr); 286 ierr = TaoLineSearchSetType(tao->linesearch, morethuente_type);CHKERRQ(ierr); 287 ierr = TaoLineSearchUseTaoRoutines(tao->linesearch,tao);CHKERRQ(ierr); 288 PetscFunctionReturn(0); 289 } 290 291