1 #include "taolinesearch.h" 2 #include "taocg.h" 3 4 #define CG_FletcherReeves 0 5 #define CG_PolakRibiere 1 6 #define CG_PolakRibierePlus 2 7 #define CG_HestenesStiefel 3 8 #define CG_DaiYuan 4 9 #define CG_Types 5 10 11 static const char *CG_Table[64] = { 12 "fr", "pr", "prp", "hs", "dy" 13 }; 14 15 16 #undef __FUNCT__ 17 #define __FUNCT__ "TaoSolve_CG" 18 static PetscErrorCode TaoSolve_CG(TaoSolver tao) 19 { 20 TAO_CG *cgP = (TAO_CG*)tao->data; 21 PetscErrorCode ierr; 22 TaoSolverTerminationReason reason = TAO_CONTINUE_ITERATING; 23 TaoLineSearchTerminationReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING; 24 PetscReal step=1.0,f,gnorm,gnorm2,delta,gd,ginner,beta; 25 PetscReal gd_old,gnorm2_old,f_old; 26 PetscInt iter=0; 27 28 29 PetscFunctionBegin; 30 31 if (tao->XL || tao->XU || tao->ops->computebounds) { 32 ierr = PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by cg algorithm\n"); CHKERRQ(ierr); 33 } 34 35 /* Check convergence criteria */ 36 ierr = TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient); CHKERRQ(ierr); 37 ierr = VecNorm(tao->gradient,NORM_2,&gnorm); CHKERRQ(ierr); 38 if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) { 39 SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN"); 40 } 41 42 ierr = TaoMonitor(tao, iter, f, gnorm, 0.0, step, &reason); CHKERRQ(ierr); 43 if (reason != TAO_CONTINUE_ITERATING) { 44 PetscFunctionReturn(0); 45 } 46 47 48 /* Set initial direction to -gradient */ 49 ierr = VecCopy(tao->gradient, tao->stepdirection); CHKERRQ(ierr); 50 ierr = VecScale(tao->stepdirection, -1.0); CHKERRQ(ierr); 51 gnorm2 = gnorm*gnorm; 52 53 /* Set initial scaling for the function */ 54 if (f != 0.0) { 55 delta = 2.0*PetscAbsScalar(f) / gnorm2; 56 delta = PetscMax(delta,cgP->delta_min); 57 delta = PetscMin(delta,cgP->delta_max); 58 } else { 59 delta = 2.0 / gnorm2; 60 delta = PetscMax(delta,cgP->delta_min); 61 delta = PetscMin(delta,cgP->delta_max); 62 } 63 /* Set counter for gradient and reset steps */ 64 cgP->ngradsteps = 0; 65 cgP->nresetsteps = 0; 66 67 while (1) { 68 /* Save the current gradient information */ 69 f_old = f; 70 gnorm2_old = gnorm2; 71 ierr = VecCopy(tao->solution, cgP->X_old); CHKERRQ(ierr); 72 ierr = VecCopy(tao->gradient, cgP->G_old); CHKERRQ(ierr); 73 ierr = VecDot(tao->gradient, tao->stepdirection, &gd); CHKERRQ(ierr); 74 if ((gd >= 0) || PetscIsInfOrNanReal(gd)) { 75 ++cgP->ngradsteps; 76 if (f != 0.0) { 77 delta = 2.0*PetscAbsScalar(f) / gnorm2; 78 delta = PetscMax(delta,cgP->delta_min); 79 delta = PetscMin(delta,cgP->delta_max); 80 } else { 81 delta = 2.0 / gnorm2; 82 delta = PetscMax(delta,cgP->delta_min); 83 delta = PetscMin(delta,cgP->delta_max); 84 } 85 86 ierr = VecCopy(tao->gradient, tao->stepdirection); CHKERRQ(ierr); 87 ierr = VecScale(tao->stepdirection, -1.0); CHKERRQ(ierr); 88 } 89 90 /* Search direction for improving point */ 91 ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,delta); 92 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status); CHKERRQ(ierr); 93 ierr = TaoAddLineSearchCounts(tao); CHKERRQ(ierr); 94 if (ls_status != TAOLINESEARCH_SUCCESS && 95 ls_status != TAOLINESEARCH_SUCCESS_USER) { 96 /* Linesearch failed */ 97 /* Reset factors and use scaled gradient step */ 98 ++cgP->nresetsteps; 99 f = f_old; 100 gnorm2 = gnorm2_old; 101 ierr = VecCopy(cgP->X_old, tao->solution); CHKERRQ(ierr); 102 ierr = VecCopy(cgP->G_old, tao->gradient); CHKERRQ(ierr); 103 104 if (f != 0.0) { 105 delta = 2.0*PetscAbsScalar(f) / gnorm2; 106 delta = PetscMax(delta,cgP->delta_min); 107 delta = PetscMin(delta,cgP->delta_max); 108 } else { 109 delta = 2.0 / gnorm2; 110 delta = PetscMax(delta,cgP->delta_min); 111 delta = PetscMin(delta,cgP->delta_max); 112 } 113 114 ierr = VecCopy(tao->gradient, tao->stepdirection); CHKERRQ(ierr); 115 ierr = VecScale(tao->stepdirection, -1.0); CHKERRQ(ierr); 116 117 ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,delta); 118 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status); CHKERRQ(ierr); 119 ierr = TaoAddLineSearchCounts(tao); CHKERRQ(ierr); 120 121 if (ls_status != TAOLINESEARCH_SUCCESS && 122 ls_status != TAOLINESEARCH_SUCCESS_USER) { 123 /* Linesearch failed again */ 124 /* switch to unscaled gradient */ 125 f = f_old; 126 gnorm2 = gnorm2_old; 127 ierr = VecCopy(cgP->X_old, tao->solution); CHKERRQ(ierr); 128 ierr = VecCopy(cgP->G_old, tao->gradient); CHKERRQ(ierr); 129 delta = 1.0; 130 ierr = VecCopy(tao->solution, tao->stepdirection); CHKERRQ(ierr); 131 ierr = VecScale(tao->stepdirection, -1.0); CHKERRQ(ierr); 132 133 ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,delta); 134 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status); CHKERRQ(ierr); 135 ierr = TaoAddLineSearchCounts(tao); CHKERRQ(ierr); 136 if (ls_status != TAOLINESEARCH_SUCCESS && 137 ls_status != TAOLINESEARCH_SUCCESS_USER) { 138 139 /* Line search failed for last time -- give up */ 140 f = f_old; 141 gnorm2 = gnorm2_old; 142 ierr = VecCopy(cgP->X_old, tao->solution); CHKERRQ(ierr); 143 ierr = VecCopy(cgP->G_old, tao->gradient); CHKERRQ(ierr); 144 step = 0.0; 145 reason = TAO_DIVERGED_LS_FAILURE; 146 tao->reason = TAO_DIVERGED_LS_FAILURE; 147 } 148 } 149 } 150 151 /* Check for bad value */ 152 ierr = VecNorm(tao->gradient,NORM_2,&gnorm); CHKERRQ(ierr); 153 if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) { 154 SETERRQ(PETSC_COMM_SELF,1,"User-provided compute function generated Inf or NaN"); 155 } 156 157 /* Check for termination */ 158 gnorm2 =gnorm * gnorm; 159 iter++; 160 ierr = TaoMonitor(tao, iter, f, gnorm, 0.0, step, &reason); CHKERRQ(ierr); 161 if (reason != TAO_CONTINUE_ITERATING) { 162 break; 163 } 164 165 /* Check for restart condition */ 166 ierr = VecDot(tao->gradient, cgP->G_old, &ginner); CHKERRQ(ierr); 167 if (PetscAbsScalar(ginner) >= cgP->eta * gnorm2) { 168 /* Gradients far from orthognal; use steepest descent direction */ 169 beta = 0.0; 170 } else { 171 /* Gradients close to orthogonal; use conjugate gradient formula */ 172 switch (cgP->cg_type) { 173 case CG_FletcherReeves: 174 beta = gnorm2 / gnorm2_old; 175 break; 176 177 case CG_PolakRibiere: 178 beta = (gnorm2 - ginner) / gnorm2_old; 179 break; 180 181 case CG_PolakRibierePlus: 182 beta = PetscMax((gnorm2-ginner)/gnorm2_old, 0.0); 183 break; 184 185 case CG_HestenesStiefel: 186 ierr = VecDot(tao->gradient, tao->stepdirection, &gd); CHKERRQ(ierr); 187 ierr = VecDot(cgP->G_old, tao->stepdirection, &gd_old); CHKERRQ(ierr); 188 beta = (gnorm2 - ginner) / (gd - gd_old); 189 break; 190 191 case CG_DaiYuan: 192 ierr = VecDot(tao->gradient, tao->stepdirection, &gd); CHKERRQ(ierr); 193 ierr = VecDot(cgP->G_old, tao->stepdirection, &gd_old); CHKERRQ(ierr); 194 beta = gnorm2 / (gd - gd_old); 195 break; 196 197 default: 198 beta = 0.0; 199 break; 200 } 201 } 202 203 /* Compute the direction d=-g + beta*d */ 204 ierr = VecAXPBY(tao->stepdirection, -1.0, beta, tao->gradient); CHKERRQ(ierr); 205 206 /* update initial steplength choice */ 207 delta = 1.0; 208 delta = PetscMax(delta, cgP->delta_min); 209 delta = PetscMin(delta, cgP->delta_max); 210 } 211 PetscFunctionReturn(0); 212 } 213 214 #undef __FUNCT__ 215 #define __FUNCT__ "TaoSetUp_CG" 216 static PetscErrorCode TaoSetUp_CG(TaoSolver tao) 217 { 218 TAO_CG *cgP = (TAO_CG*)tao->data; 219 PetscErrorCode ierr; 220 221 PetscFunctionBegin; 222 if (!tao->gradient) {ierr = VecDuplicate(tao->solution,&tao->gradient); CHKERRQ(ierr);} 223 if (!tao->stepdirection) {ierr = VecDuplicate(tao->solution,&tao->stepdirection); CHKERRQ(ierr); } 224 if (!cgP->X_old) {ierr = VecDuplicate(tao->solution,&cgP->X_old); CHKERRQ(ierr);} 225 if (!cgP->G_old) {ierr = VecDuplicate(tao->gradient,&cgP->G_old); CHKERRQ(ierr); } 226 227 PetscFunctionReturn(0); 228 } 229 230 231 #undef __FUNCT__ 232 #define __FUNCT__ "TaoDestroy_CG" 233 static PetscErrorCode TaoDestroy_CG(TaoSolver tao) 234 { 235 TAO_CG *cgP = (TAO_CG*) tao->data; 236 PetscErrorCode ierr; 237 238 PetscFunctionBegin; 239 240 if (tao->setupcalled) { 241 ierr = VecDestroy(&cgP->X_old); CHKERRQ(ierr); 242 ierr = VecDestroy(&cgP->G_old); CHKERRQ(ierr); 243 } 244 ierr = TaoLineSearchDestroy(&tao->linesearch); CHKERRQ(ierr); 245 246 ierr = PetscFree(tao->data); CHKERRQ(ierr); 247 tao->data = PETSC_NULL; 248 249 PetscFunctionReturn(0); 250 } 251 252 #undef __FUNCT__ 253 #define __FUNCT__ "TaoSetFromOptions_CG" 254 static PetscErrorCode TaoSetFromOptions_CG(TaoSolver tao) 255 { 256 TAO_CG *cgP = (TAO_CG*)tao->data; 257 PetscErrorCode ierr; 258 259 PetscFunctionBegin; 260 ierr = TaoLineSearchSetFromOptions(tao->linesearch); CHKERRQ(ierr); 261 ierr = PetscOptionsHead("Nonlinear Conjugate Gradient method for unconstrained optimization"); CHKERRQ(ierr); 262 ierr = PetscOptionsReal("-tao_cg_eta","restart tolerance", "", cgP->eta,&cgP->eta,0);CHKERRQ(ierr); 263 ierr = PetscOptionsEList("-tao_cg_type","cg formula", "", CG_Table, CG_Types, CG_Table[cgP->cg_type], &cgP->cg_type, 0); CHKERRQ(ierr); 264 ierr = PetscOptionsReal("-tao_cg_delta_min","minimum delta value", "", cgP-> 265 delta_min,&cgP->delta_min,0); CHKERRQ(ierr); 266 ierr = PetscOptionsReal("-tao_cg_delta_max","maximum delta value", "", cgP-> 267 delta_max,&cgP->delta_max,0); CHKERRQ(ierr); 268 ierr = PetscOptionsTail(); CHKERRQ(ierr); 269 PetscFunctionReturn(0); 270 } 271 272 #undef __FUNCT__ 273 #define __FUNCT__ "TaoView_CG" 274 static PetscErrorCode TaoView_CG(TaoSolver tao, PetscViewer viewer) 275 { 276 PetscBool isascii; 277 TAO_CG *cgP = (TAO_CG*)tao->data; 278 PetscErrorCode ierr; 279 280 PetscFunctionBegin; 281 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii); CHKERRQ(ierr); 282 if (isascii) { 283 ierr = PetscViewerASCIIPushTab(viewer); CHKERRQ(ierr); 284 ierr = PetscViewerASCIIPrintf(viewer, "CG Type: %s\n", CG_Table[cgP->cg_type]); CHKERRQ(ierr); 285 ierr = PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", cgP->ngradsteps); CHKERRQ(ierr); 286 ierr= PetscViewerASCIIPrintf(viewer, "Reset steps: %D\n", cgP->nresetsteps); CHKERRQ(ierr); 287 ierr = PetscViewerASCIIPopTab(viewer); CHKERRQ(ierr); 288 } else { 289 SETERRQ1(((PetscObject)tao)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for TAO CG",((PetscObject)viewer)->type_name); 290 } 291 PetscFunctionReturn(0); 292 } 293 294 295 EXTERN_C_BEGIN 296 #undef __FUNCT__ 297 #define __FUNCT__ "TaoCreate_CG" 298 PetscErrorCode TaoCreate_CG(TaoSolver tao) 299 { 300 TAO_CG *cgP; 301 const char *morethuente_type = TAOLINESEARCH_MT; 302 PetscErrorCode ierr; 303 PetscFunctionBegin; 304 tao->ops->setup = TaoSetUp_CG; 305 tao->ops->solve = TaoSolve_CG; 306 tao->ops->view = TaoView_CG; 307 tao->ops->setfromoptions = TaoSetFromOptions_CG; 308 tao->ops->destroy = TaoDestroy_CG; 309 310 tao->max_it = 2000; 311 tao->max_funcs = 4000; 312 tao->fatol = 1e-4; 313 tao->frtol = 1e-4; 314 315 /* Note: nondefault values should be used for nonlinear conjugate gradient */ 316 /* method. In particular, gtol should be less that 0.5; the value used in */ 317 /* Nocedal and Wright is 0.10. We use the default values for the */ 318 /* linesearch because it seems to work better. */ 319 ierr = TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch); CHKERRQ(ierr); 320 ierr = TaoLineSearchSetType(tao->linesearch, morethuente_type); CHKERRQ(ierr); 321 ierr = TaoLineSearchUseTaoSolverRoutines(tao->linesearch, tao); CHKERRQ(ierr); 322 323 324 ierr = PetscNewLog(tao,&cgP); CHKERRQ(ierr); 325 tao->data = (void*)cgP; 326 cgP->eta = 0.1; 327 cgP->delta_min = 1e-7; 328 cgP->delta_max = 100; 329 cgP->cg_type = CG_PolakRibierePlus; 330 331 PetscFunctionReturn(0); 332 } 333 334 EXTERN_C_END 335 336 337 338 339 340