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