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 if (tao->setupcalled) { 240 ierr = VecDestroy(&cgP->X_old); CHKERRQ(ierr); 241 ierr = VecDestroy(&cgP->G_old); CHKERRQ(ierr); 242 } 243 ierr = TaoLineSearchDestroy(&tao->linesearch); CHKERRQ(ierr); 244 ierr = PetscFree(tao->data); CHKERRQ(ierr); 245 PetscFunctionReturn(0); 246 } 247 248 #undef __FUNCT__ 249 #define __FUNCT__ "TaoSetFromOptions_CG" 250 static PetscErrorCode TaoSetFromOptions_CG(TaoSolver tao) 251 { 252 TAO_CG *cgP = (TAO_CG*)tao->data; 253 PetscErrorCode ierr; 254 255 PetscFunctionBegin; 256 ierr = TaoLineSearchSetFromOptions(tao->linesearch); CHKERRQ(ierr); 257 ierr = PetscOptionsHead("Nonlinear Conjugate Gradient method for unconstrained optimization"); CHKERRQ(ierr); 258 ierr = PetscOptionsReal("-tao_cg_eta","restart tolerance", "", cgP->eta,&cgP->eta,0);CHKERRQ(ierr); 259 ierr = PetscOptionsEList("-tao_cg_type","cg formula", "", CG_Table, CG_Types, CG_Table[cgP->cg_type], &cgP->cg_type, 0); CHKERRQ(ierr); 260 ierr = PetscOptionsReal("-tao_cg_delta_min","minimum delta value", "", cgP->delta_min,&cgP->delta_min,0); CHKERRQ(ierr); 261 ierr = PetscOptionsReal("-tao_cg_delta_max","maximum delta value", "", cgP->delta_max,&cgP->delta_max,0); CHKERRQ(ierr); 262 ierr = PetscOptionsTail(); CHKERRQ(ierr); 263 PetscFunctionReturn(0); 264 } 265 266 #undef __FUNCT__ 267 #define __FUNCT__ "TaoView_CG" 268 static PetscErrorCode TaoView_CG(TaoSolver tao, PetscViewer viewer) 269 { 270 PetscBool isascii; 271 TAO_CG *cgP = (TAO_CG*)tao->data; 272 PetscErrorCode ierr; 273 274 PetscFunctionBegin; 275 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii); CHKERRQ(ierr); 276 if (isascii) { 277 ierr = PetscViewerASCIIPushTab(viewer); CHKERRQ(ierr); 278 ierr = PetscViewerASCIIPrintf(viewer, "CG Type: %s\n", CG_Table[cgP->cg_type]); CHKERRQ(ierr); 279 ierr = PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", cgP->ngradsteps); CHKERRQ(ierr); 280 ierr= PetscViewerASCIIPrintf(viewer, "Reset steps: %D\n", cgP->nresetsteps); CHKERRQ(ierr); 281 ierr = PetscViewerASCIIPopTab(viewer); CHKERRQ(ierr); 282 } else { 283 SETERRQ1(((PetscObject)tao)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for TAO CG",((PetscObject)viewer)->type_name); 284 } 285 PetscFunctionReturn(0); 286 } 287 288 289 EXTERN_C_BEGIN 290 #undef __FUNCT__ 291 #define __FUNCT__ "TaoCreate_CG" 292 PetscErrorCode TaoCreate_CG(TaoSolver tao) 293 { 294 TAO_CG *cgP; 295 const char *morethuente_type = TAOLINESEARCH_MT; 296 PetscErrorCode ierr; 297 PetscFunctionBegin; 298 tao->ops->setup = TaoSetUp_CG; 299 tao->ops->solve = TaoSolve_CG; 300 tao->ops->view = TaoView_CG; 301 tao->ops->setfromoptions = TaoSetFromOptions_CG; 302 tao->ops->destroy = TaoDestroy_CG; 303 304 tao->max_it = 2000; 305 tao->max_funcs = 4000; 306 tao->fatol = 1e-4; 307 tao->frtol = 1e-4; 308 309 /* Note: nondefault values should be used for nonlinear conjugate gradient */ 310 /* method. In particular, gtol should be less that 0.5; the value used in */ 311 /* Nocedal and Wright is 0.10. We use the default values for the */ 312 /* linesearch because it seems to work better. */ 313 ierr = TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch); CHKERRQ(ierr); 314 ierr = TaoLineSearchSetType(tao->linesearch, morethuente_type); CHKERRQ(ierr); 315 ierr = TaoLineSearchUseTaoSolverRoutines(tao->linesearch, tao); CHKERRQ(ierr); 316 317 318 ierr = PetscNewLog(tao,&cgP); CHKERRQ(ierr); 319 tao->data = (void*)cgP; 320 cgP->eta = 0.1; 321 cgP->delta_min = 1e-7; 322 cgP->delta_max = 100; 323 cgP->cg_type = CG_PolakRibierePlus; 324 325 PetscFunctionReturn(0); 326 } 327 328 EXTERN_C_END 329 330 331 332 333 334