1 #include <petscksp.h> 2 #include <../src/tao/quadratic/impls/gpcg/gpcg.h> /*I "gpcg.h" I*/ 3 4 static PetscErrorCode GPCGGradProjections(Tao tao); 5 static PetscErrorCode GPCGObjectiveAndGradient(TaoLineSearch,Vec,PetscReal*,Vec,void*); 6 7 /*------------------------------------------------------------*/ 8 static PetscErrorCode TaoDestroy_GPCG(Tao tao) 9 { 10 TAO_GPCG *gpcg = (TAO_GPCG *)tao->data; 11 12 /* Free allocated memory in GPCG structure */ 13 PetscFunctionBegin; 14 PetscCall(VecDestroy(&gpcg->B)); 15 PetscCall(VecDestroy(&gpcg->Work)); 16 PetscCall(VecDestroy(&gpcg->X_New)); 17 PetscCall(VecDestroy(&gpcg->G_New)); 18 PetscCall(VecDestroy(&gpcg->DXFree)); 19 PetscCall(VecDestroy(&gpcg->R)); 20 PetscCall(VecDestroy(&gpcg->PG)); 21 PetscCall(MatDestroy(&gpcg->Hsub)); 22 PetscCall(MatDestroy(&gpcg->Hsub_pre)); 23 PetscCall(ISDestroy(&gpcg->Free_Local)); 24 PetscCall(PetscFree(tao->data)); 25 PetscFunctionReturn(0); 26 } 27 28 /*------------------------------------------------------------*/ 29 static PetscErrorCode TaoSetFromOptions_GPCG(PetscOptionItems *PetscOptionsObject,Tao tao) 30 { 31 TAO_GPCG *gpcg = (TAO_GPCG *)tao->data; 32 PetscBool flg; 33 34 PetscFunctionBegin; 35 PetscOptionsHeadBegin(PetscOptionsObject,"Gradient Projection, Conjugate Gradient method for bound constrained optimization"); 36 PetscCall(PetscOptionsInt("-tao_gpcg_maxpgits","maximum number of gradient projections per GPCG iterate",NULL,gpcg->maxgpits,&gpcg->maxgpits,&flg)); 37 PetscOptionsHeadEnd(); 38 PetscCall(KSPSetFromOptions(tao->ksp)); 39 PetscCall(TaoLineSearchSetFromOptions(tao->linesearch)); 40 PetscFunctionReturn(0); 41 } 42 43 /*------------------------------------------------------------*/ 44 static PetscErrorCode TaoView_GPCG(Tao tao, PetscViewer viewer) 45 { 46 TAO_GPCG *gpcg = (TAO_GPCG *)tao->data; 47 PetscBool isascii; 48 49 PetscFunctionBegin; 50 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii)); 51 if (isascii) { 52 PetscCall(PetscViewerASCIIPrintf(viewer,"Total PG its: %" PetscInt_FMT ",",gpcg->total_gp_its)); 53 PetscCall(PetscViewerASCIIPrintf(viewer,"PG tolerance: %g \n",(double)gpcg->pg_ftol)); 54 } 55 PetscCall(TaoLineSearchView(tao->linesearch,viewer)); 56 PetscFunctionReturn(0); 57 } 58 59 /* GPCGObjectiveAndGradient() 60 Compute f=0.5 * x'Hx + b'x + c 61 g=Hx + b 62 */ 63 static PetscErrorCode GPCGObjectiveAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void*tptr) 64 { 65 Tao tao = (Tao)tptr; 66 TAO_GPCG *gpcg = (TAO_GPCG*)tao->data; 67 PetscReal f1,f2; 68 69 PetscFunctionBegin; 70 PetscCall(MatMult(tao->hessian,X,G)); 71 PetscCall(VecDot(G,X,&f1)); 72 PetscCall(VecDot(gpcg->B,X,&f2)); 73 PetscCall(VecAXPY(G,1.0,gpcg->B)); 74 *f=f1/2.0 + f2 + gpcg->c; 75 PetscFunctionReturn(0); 76 } 77 78 /* ---------------------------------------------------------- */ 79 static PetscErrorCode TaoSetup_GPCG(Tao tao) 80 { 81 TAO_GPCG *gpcg = (TAO_GPCG *)tao->data; 82 83 PetscFunctionBegin; 84 /* Allocate some arrays */ 85 if (!tao->gradient) { 86 PetscCall(VecDuplicate(tao->solution, &tao->gradient)); 87 } 88 if (!tao->stepdirection) { 89 PetscCall(VecDuplicate(tao->solution, &tao->stepdirection)); 90 } 91 if (!tao->XL) { 92 PetscCall(VecDuplicate(tao->solution,&tao->XL)); 93 PetscCall(VecSet(tao->XL,PETSC_NINFINITY)); 94 } 95 if (!tao->XU) { 96 PetscCall(VecDuplicate(tao->solution,&tao->XU)); 97 PetscCall(VecSet(tao->XU,PETSC_INFINITY)); 98 } 99 100 PetscCall(VecDuplicate(tao->solution,&gpcg->B)); 101 PetscCall(VecDuplicate(tao->solution,&gpcg->Work)); 102 PetscCall(VecDuplicate(tao->solution,&gpcg->X_New)); 103 PetscCall(VecDuplicate(tao->solution,&gpcg->G_New)); 104 PetscCall(VecDuplicate(tao->solution,&gpcg->DXFree)); 105 PetscCall(VecDuplicate(tao->solution,&gpcg->R)); 106 PetscCall(VecDuplicate(tao->solution,&gpcg->PG)); 107 /* 108 if (gpcg->ksp_type == GPCG_KSP_NASH) { 109 PetscCall(KSPSetType(tao->ksp,KSPNASH)); 110 } else if (gpcg->ksp_type == GPCG_KSP_STCG) { 111 PetscCall(KSPSetType(tao->ksp,KSPSTCG)); 112 } else { 113 PetscCall(KSPSetType(tao->ksp,KSPGLTR)); 114 } 115 if (tao->ksp->ops->setfromoptions) { 116 (*tao->ksp->ops->setfromoptions)(tao->ksp); 117 } 118 119 } 120 */ 121 PetscFunctionReturn(0); 122 } 123 124 static PetscErrorCode TaoSolve_GPCG(Tao tao) 125 { 126 TAO_GPCG *gpcg = (TAO_GPCG *)tao->data; 127 PetscInt its; 128 PetscReal actred,f,f_new,gnorm,gdx,stepsize,xtb; 129 PetscReal xtHx; 130 TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING; 131 132 PetscFunctionBegin; 133 134 PetscCall(TaoComputeVariableBounds(tao)); 135 PetscCall(VecMedian(tao->XL,tao->solution,tao->XU,tao->solution)); 136 PetscCall(TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU)); 137 138 /* Using f = .5*x'Hx + x'b + c and g=Hx + b, compute b,c */ 139 PetscCall(TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre)); 140 PetscCall(TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient)); 141 PetscCall(VecCopy(tao->gradient, gpcg->B)); 142 PetscCall(MatMult(tao->hessian,tao->solution,gpcg->Work)); 143 PetscCall(VecDot(gpcg->Work, tao->solution, &xtHx)); 144 PetscCall(VecAXPY(gpcg->B,-1.0,gpcg->Work)); 145 PetscCall(VecDot(gpcg->B,tao->solution,&xtb)); 146 gpcg->c=f-xtHx/2.0-xtb; 147 if (gpcg->Free_Local) { 148 PetscCall(ISDestroy(&gpcg->Free_Local)); 149 } 150 PetscCall(VecWhichInactive(tao->XL,tao->solution,tao->gradient,tao->XU,PETSC_TRUE,&gpcg->Free_Local)); 151 152 /* Project the gradient and calculate the norm */ 153 PetscCall(VecCopy(tao->gradient,gpcg->G_New)); 154 PetscCall(VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,gpcg->PG)); 155 PetscCall(VecNorm(gpcg->PG,NORM_2,&gpcg->gnorm)); 156 tao->step=1.0; 157 gpcg->f = f; 158 159 /* Check Stopping Condition */ 160 tao->reason = TAO_CONTINUE_ITERATING; 161 PetscCall(TaoLogConvergenceHistory(tao,f,gpcg->gnorm,0.0,tao->ksp_its)); 162 PetscCall(TaoMonitor(tao,tao->niter,f,gpcg->gnorm,0.0,tao->step)); 163 PetscCall((*tao->ops->convergencetest)(tao,tao->cnvP)); 164 165 while (tao->reason == TAO_CONTINUE_ITERATING) { 166 /* Call general purpose update function */ 167 if (tao->ops->update) { 168 PetscCall((*tao->ops->update)(tao, tao->niter, tao->user_update)); 169 } 170 tao->ksp_its=0; 171 172 PetscCall(GPCGGradProjections(tao)); 173 PetscCall(ISGetSize(gpcg->Free_Local,&gpcg->n_free)); 174 175 f=gpcg->f; gnorm=gpcg->gnorm; 176 177 PetscCall(KSPReset(tao->ksp)); 178 179 if (gpcg->n_free > 0) { 180 /* Create a reduced linear system */ 181 PetscCall(VecDestroy(&gpcg->R)); 182 PetscCall(VecDestroy(&gpcg->DXFree)); 183 PetscCall(TaoVecGetSubVec(tao->gradient,gpcg->Free_Local, tao->subset_type, 0.0, &gpcg->R)); 184 PetscCall(VecScale(gpcg->R, -1.0)); 185 PetscCall(TaoVecGetSubVec(tao->stepdirection,gpcg->Free_Local,tao->subset_type, 0.0, &gpcg->DXFree)); 186 PetscCall(VecSet(gpcg->DXFree,0.0)); 187 188 PetscCall(TaoMatGetSubMat(tao->hessian, gpcg->Free_Local, gpcg->Work, tao->subset_type, &gpcg->Hsub)); 189 190 if (tao->hessian_pre == tao->hessian) { 191 PetscCall(MatDestroy(&gpcg->Hsub_pre)); 192 PetscCall(PetscObjectReference((PetscObject)gpcg->Hsub)); 193 gpcg->Hsub_pre = gpcg->Hsub; 194 } else { 195 PetscCall(TaoMatGetSubMat(tao->hessian, gpcg->Free_Local, gpcg->Work, tao->subset_type, &gpcg->Hsub_pre)); 196 } 197 198 PetscCall(KSPReset(tao->ksp)); 199 PetscCall(KSPSetOperators(tao->ksp,gpcg->Hsub,gpcg->Hsub_pre)); 200 201 PetscCall(KSPSolve(tao->ksp,gpcg->R,gpcg->DXFree)); 202 PetscCall(KSPGetIterationNumber(tao->ksp,&its)); 203 tao->ksp_its+=its; 204 tao->ksp_tot_its+=its; 205 PetscCall(VecSet(tao->stepdirection,0.0)); 206 PetscCall(VecISAXPY(tao->stepdirection,gpcg->Free_Local,1.0,gpcg->DXFree)); 207 208 PetscCall(VecDot(tao->stepdirection,tao->gradient,&gdx)); 209 PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch,1.0)); 210 f_new=f; 211 PetscCall(TaoLineSearchApply(tao->linesearch,tao->solution,&f_new,tao->gradient,tao->stepdirection,&stepsize,&ls_status)); 212 213 actred = f_new - f; 214 215 /* Evaluate the function and gradient at the new point */ 216 PetscCall(VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU, gpcg->PG)); 217 PetscCall(VecNorm(gpcg->PG, NORM_2, &gnorm)); 218 f=f_new; 219 PetscCall(ISDestroy(&gpcg->Free_Local)); 220 PetscCall(VecWhichInactive(tao->XL,tao->solution,tao->gradient,tao->XU,PETSC_TRUE,&gpcg->Free_Local)); 221 } else { 222 actred = 0; gpcg->step=1.0; 223 /* if there were no free variables, no cg method */ 224 } 225 226 tao->niter++; 227 gpcg->f=f;gpcg->gnorm=gnorm; gpcg->actred=actred; 228 PetscCall(TaoLogConvergenceHistory(tao,f,gpcg->gnorm,0.0,tao->ksp_its)); 229 PetscCall(TaoMonitor(tao,tao->niter,f,gpcg->gnorm,0.0,tao->step)); 230 PetscCall((*tao->ops->convergencetest)(tao,tao->cnvP)); 231 if (tao->reason != TAO_CONTINUE_ITERATING) break; 232 } /* END MAIN LOOP */ 233 234 PetscFunctionReturn(0); 235 } 236 237 static PetscErrorCode GPCGGradProjections(Tao tao) 238 { 239 TAO_GPCG *gpcg = (TAO_GPCG *)tao->data; 240 PetscInt i; 241 PetscReal actred=-1.0,actred_max=0.0, gAg,gtg=gpcg->gnorm,alpha; 242 PetscReal f_new,gdx,stepsize; 243 Vec DX=tao->stepdirection,XL=tao->XL,XU=tao->XU,Work=gpcg->Work; 244 Vec X=tao->solution,G=tao->gradient; 245 TaoLineSearchConvergedReason lsflag=TAOLINESEARCH_CONTINUE_ITERATING; 246 247 /* 248 The free, active, and binding variables should be already identified 249 */ 250 PetscFunctionBegin; 251 for (i=0;i<gpcg->maxgpits;i++) { 252 if (-actred <= (gpcg->pg_ftol)*actred_max) break; 253 PetscCall(VecBoundGradientProjection(G,X,XL,XU,DX)); 254 PetscCall(VecScale(DX,-1.0)); 255 PetscCall(VecDot(DX,G,&gdx)); 256 257 PetscCall(MatMult(tao->hessian,DX,Work)); 258 PetscCall(VecDot(DX,Work,&gAg)); 259 260 gpcg->gp_iterates++; 261 gpcg->total_gp_its++; 262 263 gtg=-gdx; 264 if (PetscAbsReal(gAg) == 0.0) { 265 alpha = 1.0; 266 } else { 267 alpha = PetscAbsReal(gtg/gAg); 268 } 269 PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch,alpha)); 270 f_new=gpcg->f; 271 PetscCall(TaoLineSearchApply(tao->linesearch,X,&f_new,G,DX,&stepsize,&lsflag)); 272 273 /* Update the iterate */ 274 actred = f_new - gpcg->f; 275 actred_max = PetscMax(actred_max,-(f_new - gpcg->f)); 276 gpcg->f = f_new; 277 PetscCall(ISDestroy(&gpcg->Free_Local)); 278 PetscCall(VecWhichInactive(XL,X,tao->gradient,XU,PETSC_TRUE,&gpcg->Free_Local)); 279 } 280 281 gpcg->gnorm=gtg; 282 PetscFunctionReturn(0); 283 } /* End gradient projections */ 284 285 static PetscErrorCode TaoComputeDual_GPCG(Tao tao, Vec DXL, Vec DXU) 286 { 287 TAO_GPCG *gpcg = (TAO_GPCG *)tao->data; 288 289 PetscFunctionBegin; 290 PetscCall(VecBoundGradientProjection(tao->gradient, tao->solution, tao->XL, tao->XU, gpcg->Work)); 291 PetscCall(VecCopy(gpcg->Work, DXL)); 292 PetscCall(VecAXPY(DXL,-1.0,tao->gradient)); 293 PetscCall(VecSet(DXU,0.0)); 294 PetscCall(VecPointwiseMax(DXL,DXL,DXU)); 295 296 PetscCall(VecCopy(tao->gradient,DXU)); 297 PetscCall(VecAXPY(DXU,-1.0,gpcg->Work)); 298 PetscCall(VecSet(gpcg->Work,0.0)); 299 PetscCall(VecPointwiseMin(DXU,gpcg->Work,DXU)); 300 PetscFunctionReturn(0); 301 } 302 303 /*------------------------------------------------------------*/ 304 /*MC 305 TAOGPCG - gradient projected conjugate gradient algorithm is an active-set 306 conjugate-gradient based method for bound-constrained minimization 307 308 Options Database Keys: 309 + -tao_gpcg_maxpgits - maximum number of gradient projections for GPCG iterate 310 - -tao_subset_type - "subvec","mask","matrix-free", strategies for handling active-sets 311 312 Level: beginner 313 M*/ 314 PETSC_EXTERN PetscErrorCode TaoCreate_GPCG(Tao tao) 315 { 316 TAO_GPCG *gpcg; 317 318 PetscFunctionBegin; 319 tao->ops->setup = TaoSetup_GPCG; 320 tao->ops->solve = TaoSolve_GPCG; 321 tao->ops->view = TaoView_GPCG; 322 tao->ops->setfromoptions = TaoSetFromOptions_GPCG; 323 tao->ops->destroy = TaoDestroy_GPCG; 324 tao->ops->computedual = TaoComputeDual_GPCG; 325 326 PetscCall(PetscNewLog(tao,&gpcg)); 327 tao->data = (void*)gpcg; 328 329 /* Override default settings (unless already changed) */ 330 if (!tao->max_it_changed) tao->max_it=500; 331 if (!tao->max_funcs_changed) tao->max_funcs = 100000; 332 #if defined(PETSC_USE_REAL_SINGLE) 333 if (!tao->gatol_changed) tao->gatol=1e-6; 334 if (!tao->grtol_changed) tao->grtol=1e-6; 335 #else 336 if (!tao->gatol_changed) tao->gatol=1e-12; 337 if (!tao->grtol_changed) tao->grtol=1e-12; 338 #endif 339 340 /* Initialize pointers and variables */ 341 gpcg->n=0; 342 gpcg->maxgpits = 8; 343 gpcg->pg_ftol = 0.1; 344 345 gpcg->gp_iterates=0; /* Cumulative number */ 346 gpcg->total_gp_its = 0; 347 348 /* Initialize pointers and variables */ 349 gpcg->n_bind=0; 350 gpcg->n_free = 0; 351 gpcg->n_upper=0; 352 gpcg->n_lower=0; 353 gpcg->subset_type = TAO_SUBSET_MASK; 354 gpcg->Hsub=NULL; 355 gpcg->Hsub_pre=NULL; 356 357 PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp)); 358 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1)); 359 PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix)); 360 PetscCall(KSPSetType(tao->ksp,KSPNASH)); 361 362 PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch)); 363 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1)); 364 PetscCall(TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHGPCG)); 365 PetscCall(TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch, GPCGObjectiveAndGradient, tao)); 366 PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix)); 367 PetscFunctionReturn(0); 368 } 369