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