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 /*------------------------------------------------------------*/
TaoDestroy_GPCG(Tao tao)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 /*------------------------------------------------------------*/
TaoSetFromOptions_GPCG(PetscOptionItems * PetscOptionsObject,Tao tao)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 /*------------------------------------------------------------*/
TaoView_GPCG(Tao tao,PetscViewer viewer)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 */
GPCGObjectiveAndGradient(TaoLineSearch ls,Vec X,PetscReal * f,Vec G,void * tptr)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 /* ---------------------------------------------------------- */
TaoSetup_GPCG(Tao tao)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,KSPNASH);CHKERRQ(ierr);
116       } else if (gpcg->ksp_type == GPCG_KSP_STCG) {
117         ierr = KSPSetType(tao->ksp,KSPSTCG);CHKERRQ(ierr);
118       } else {
119         ierr = KSPSetType(tao->ksp,KSPGLTR);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 
TaoSolve_GPCG(Tao tao)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 
GPCGGradProjections(Tao tao)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 
TaoComputeDual_GPCG(Tao tao,Vec DXL,Vec DXU)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*/
TaoCreate_GPCG(Tao tao)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,KSPNASH);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