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