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