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