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