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