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