Actual source code: daapp_mgrid.c

  1: #include "taodaapplication.h"     /* taodaapplication.h */
  2: //#include "taoapp.h"
 3:  #include src/petsctao/include/taopetsc.h

  5: #include "../../application/taoapp/taoapp_petsc.h"
  6: #include "../interface/daapp_impl.h"
 7:  #include multigridmat.h

  9: #include "petscmg.h"

 11: class TaoMultiGridApplication: public TaoPetscApplication {

 13:  protected:

 15:  public:

 17:   TaoMultiGridApplication(MPI_Comm);
 18:   TaoMultiGridApplication();

 20:   int coarselevel;
 21:   int EvaluateHessian(TaoVec *xx, TaoMat *HH);

 23: };


 28: TaoMultiGridApplication :: TaoMultiGridApplication(MPI_Comm mpicomm) : TaoPetscApplication(mpicomm) {
 29:   Mat M=0;
 30:   this->coarselevel=0;
 31:   this->taoh=new TaoMultiGridMatPetsc(M);
 32:   this->taoj=new TaoMultiGridMatPetsc(M);
 33:   return;
 34: }

 38: TaoMultiGridApplication::TaoMultiGridApplication() : TaoPetscApplication(PETSC_COMM_WORLD){
 39: 
 40:   Mat M=0;
 41:   this->coarselevel=0;
 42:   this->taoh=new TaoMultiGridMatPetsc(M);
 43:   this->taoj=new TaoMultiGridMatPetsc(M);
 44:   return;
 45: }


 50: int TaoMultiGridApplication::EvaluateHessian(TaoVec *xx, TaoMat *HH){
 51:   int     info,i,clevel,cclevel;
 52:   Vec X = (Vec) xx->VecObject;
 53:   TaoMatPetsc *MM=(TaoMatPetsc*)HH;
 54:   Mat H, HPre;
 55:   MatStructure pflag;
 56:   DA_APPLICATION daapp;

 59:   info=TaoAppDAApp(this->papp,&daapp); CHKERRQ(info);

 61:   PetscStackPush("TAO User DA Hessian");
 62:   info=MM->GetMatrix(&H,&HPre,&pflag);CHKERRQ(info);
 63:   info=DAAppGetCurrentLevel(this->papp,&clevel); CHKERRQ(info);
 64:   cclevel=PetscMin(clevel,this->coarselevel);
 65:   cclevel=PetscMax(cclevel,0);

 67:   for (i=clevel; i>=cclevel; i--){
 68:     daapp->currentlevel=i;
 69:     info=TaoAppSetColoring(this->papp,daapp->grid[i].coloring); CHKERRQ(info);
 70:     info = PetscLogInfo((daapp,"TAO hessian evaluation at DA_APPLICATION object, level %d.\n",i)); CHKERRQ(info);
 71:     info = TaoAppComputeHessian(this->papp, X, &H, &HPre, &pflag); CHKERRQ(info);
 72:     if (i==clevel){
 73:       info = MM->SetMatrix(H,HPre,pflag);CHKERRQ(info);
 74:     }

 76:     if (i>cclevel){
 77:       info=MatMultTranspose(daapp->grid[i].Interpolate,X,daapp->grid[i-1].R);CHKERRQ(info);
 78:       info=VecPointwiseMult(daapp->grid[i-1].R,daapp->grid[i].CScale,daapp->grid[i-1].R);CHKERRQ(info);
 79:       X=daapp->grid[i-1].R;
 80:       H=daapp->grid[i-1].H;
 81:       HPre=H;
 82:     }
 83:   }
 84:   daapp->currentlevel=clevel;
 85:   info=TaoAppSetColoring(this->papp,daapp->grid[clevel].coloring); CHKERRQ(info);
 86:   PetscStackPop;

 88:   return(0);
 89: }

 91: int DAAppSetMultigridKSP(GridCtx *, int, KSP);
 92: int DAAppSetupMultigridMonitor(TAO_APPLICATION, DA, int, void *);
 93: int TaoAppGetMultiGridApplication(TAO_APPLICATION, TaoMultiGridApplication *);
 94: static char TaoPetscPCMGDAAppXXX[] = "TaoPetscPCMGDAApp";
 95: static char TaoPetscAppXXX[] = "TaoPetscApp";


100: static int DAAppDestroyTaoAppXXX(void*ctx){
101:   TaoMultiGridApplication *mctx=(TaoMultiGridApplication*)ctx;
103:   if (mctx){ delete mctx;}
104:   return(0);
105: }

109: int TaoAppGetMultiGridApplication(TAO_APPLICATION daapplication, TaoMultiGridApplication **mgdaapp){
110:   int ii,clevel,info;
111:   Mat H,HH;
112:   MatStructure         pflag=SAME_NONZERO_PATTERN;
113:   DA_APPLICATION       daapp;
114:   TaoMultiGridApplication*   daapp2;
115:   TaoPetscApplication  *pppappp;

119:   info=TaoAppQueryForObject(daapplication,TaoPetscPCMGDAAppXXX,(void**)&daapp2); CHKERRQ(info);
120:   if (daapp2==0){
121:     daapp2=new TaoMultiGridApplication(PETSC_COMM_WORLD);
122:     daapp2->papp=daapplication;
123:     info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);

125:     info=TaoAppQueryForObject(daapplication,TaoPetscAppXXX,(void**)&pppappp); CHKERRQ(info);
126:     if (pppappp && pppappp->taoh && 0==1){
127:       info=pppappp->taoh->GetMatrix(&H,&HH,&pflag); CHKERRQ(info);
128:       info=daapp2->taoh->SetMatrix(H,HH,pflag);CHKERRQ(info);
129:     } else {
130:       info=DAAppGetCurrentLevel(daapplication,&clevel); CHKERRQ(info);
131:       info=daapp2->taoh->SetMatrix(daapp->grid[clevel].H,daapp->grid[clevel].H,pflag);CHKERRQ(info);
132:     }
133:     if (pppappp){
134:       daapp2->tao=pppappp->tao;
135:     }
136:     info=TaoAppAddObject(daapplication,TaoPetscPCMGDAAppXXX,(void*)daapp2,&ii); CHKERRQ(info);
137:     info=TaoAppQueryRemoveObject(daapplication,TaoPetscAppXXX); CHKERRQ(info);
138:     info=TaoAppAddObject(daapplication,TaoPetscAppXXX,(void*)daapp2,&ii); CHKERRQ(info);
139:     info=TaoAppSetDestroyRoutine(daapplication,DAAppDestroyTaoAppXXX, (void*)daapp2); CHKERRQ(info);
140:   }
141:   *mgdaapp=daapp2;
142:   return(0);
143: }

145: static int UsePCMGGG=0;
148: /*@
149:   DAAppUseMultigrid - Set the preconditioner for the linear solver to be an algebraic multigrid.

151:   Collective on TAO_APPLICATION
152:   
153:   Input Parameters:
154: +  daapplication - the DA Application object
155: -  coarselevel - the coarsest grid to be used in the multigrid preconditioner. (Grid 0 is the coarsest grid.

157:    Level: intermediate

159:    Options Database Key:
160: +  -tao_da_multigrid - use multigrid linear solver
161: -  -ksp_view - view the linear solver

163: Note:
164:   This function should be called after DAAppSetHessianRoutine();

166: Note:
167:   This function should be called before each optimization solver as part of the DAAppMonitor

169: Note:
170:   Multigrid functionality is still under developement for good performance.

172: .seealso: TaoGetKSP(), DAAppSetupMultigrid()

174:    Options Database Key:
175: .  -tao_da_multigrid

177: .keywords: Linear Solver, Multigrid, DA, KSP

179: @*/
180: int DAAppUseMultigrid(TAO_APPLICATION daapplication, int coarselevel){
181:   int info;
183:   UsePCMGGG=coarselevel;
184:   info=DAAppSetBeforeMonitor(daapplication,DAAppSetupMultigridMonitor,0);
185:   CHKERRQ(info);
186:   return(0);
187: }


192: int DAAppSetupMultigridMonitor(TAO_APPLICATION daapplication, DA da, int level, void *ctx){
193:   int info;
195:   info = DAAppSetupMultigrid(daapplication,UsePCMGGG); CHKERRQ(info);
196:   return(0);
197: }


202: /*@
203:    DAAppSetupMultigrid - Sets up the multigrid preconditioner.

205:    Collective on TAO_APPLICATION

207:    Input Parameters:
208: +  daapplication - the TAO_APPLICATION context
209: -  coarselevel - the coarsest grid that should be used in the multigrid (>=0)

211:    Note:
212:    Usually the coarselevel is set to 0;

214:    Level: intermediate

216: .seealso:  DAAppUseMultigrid(), TaoAppSetMonitor()
217: @*/
218: int DAAppSetupMultigrid(TAO_APPLICATION daapplication, int coarselevel){
219:   int nn,level,info;
220:   TaoMultiGridMatPetsc *MMM;
221:   DA_APPLICATION daapp;
222:   TaoMultiGridApplication *mgdaapp;
223:   TAO_SOLVER tao;
224:   KSP  ksp;
225: 
227:   info=DAAppGetCurrentLevel(daapplication,&level); CHKERRQ(info);
228:   if (coarselevel<level){
229:     info=TaoAppGetMultiGridApplication(daapplication,&mgdaapp); CHKERRQ(info);
230:     info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);
231:     mgdaapp->coarselevel=PetscMax(0,coarselevel);
232:     nn=level-coarselevel+1;
233:     if ( mgdaapp->taoh ){
234:       MMM=(TaoMultiGridMatPetsc*)(mgdaapp->taoh);
235:       info=MMM->SetUpMultiGrid(daapp->grid+coarselevel,nn); CHKERRQ(info);
236:     } else if ( mgdaapp->taoj ){
237:       MMM=(TaoMultiGridMatPetsc*)(mgdaapp->taoj);
238:       info=MMM->SetUpMultiGrid(daapp->grid+coarselevel,nn); CHKERRQ(info);
239:       mgdaapp->taoj=MMM;
240:     }

242:     info=TaoAppGetTaoSolver(daapplication,&tao); CHKERRQ(info);
243:     info=TaoSetupApplicationSolver(daapplication,tao); CHKERRQ(info);
244:     info=TaoAppGetKSP(daapplication,&ksp); CHKERRQ(info);
245:     if (ksp){
247:       info = DAAppSetMultigridKSP(daapp->grid+coarselevel,nn,ksp);CHKERRQ(info);
248:     }
249:   }
250:   return(0);
251: }


254: 
257: int DAAppSetMultigridKSP(GridCtx* dagrid, int nlevels, KSP ksp){
258:   int i,info;
259:   PC pc;
260:   KSP sksp;
261:   KSPType ksptype;

264:   if (ksp==0 || nlevels<=1){
265:     return(0);
266:   }

269:   info = PetscLogInfo((ksp,"Set multigrid precondition in optimization solver.\n")); CHKERRQ(info);
270: 
271:   info = KSPGetType(ksp,&ksptype); CHKERRQ(info);

273:   info = KSPGetPC(ksp,&pc); CHKERRQ(info);
274:   info = PCSetType(pc,PCMG); CHKERRQ(info);
275: 
276:   info = PCMGSetLevels(pc,nlevels,PETSC_NULL); CHKERRQ(info);
277: 
278:   for (i=0; i<nlevels; i++) {
279: 
280:     info = PCMGGetSmoother(pc,i,&sksp); CHKERRQ(info);
281:     info = KSPSetType(sksp,ksptype); CHKERRQ(info);

283:     info = KSPSetOperators(sksp,dagrid[i].H,dagrid[i].H,SAME_NONZERO_PATTERN); CHKERRQ(info);
284: 
285:     info = PCMGSetX(pc,i,dagrid[i].X); CHKERRQ(info);
286:     info = PCMGSetRhs(pc,i,dagrid[i].RHS); CHKERRQ(info);
287:     info = PCMGSetR(pc,i,dagrid[i].R); CHKERRQ(info);
288:     info = PCMGSetResidual(pc,i,PCMGDefaultResidual,dagrid[i].H); CHKERRQ(info);
289:   }
290:   for (i=1; i<nlevels; i++) {
291:     info = PCMGSetInterpolate(pc,i,dagrid[i].Interpolate); CHKERRQ(info);
292:     info = PCMGSetRestriction(pc,i,dagrid[i].Interpolate); CHKERRQ(info);
293:   }
294: 
295:   return(0);
296: }

300: /* @
301:   DAAppSetQuadraticObjective - identify the objective function as quadratic or not.

303:   Collective on TAO_APPLICATION

305:   Input Parameters:
306: +  daapplication - the DA Application object
307: -  flag - indicates whether the objective is quadratic or not

309:    Level: intermediate

311:    Options Database Key:
312: .  -tao_da_quadratic - use multigrid linear solver

314:    Note:
315:    The default value is PETSC_FALSE.

317:    Note:
318:    If quadratic, consider setting the flag used in KSP to SAME_PRECONDITIONER

320:    Note:
321:    If quadratic, this routine reduces the number of Hessian evaluations done when using
322:    the multigrid preconditioner.

324: .seealso: DAAppSetMatStructure(),  DAAppUseMultigrid()


327: .keywords: DA, Objective Function

329: @ */
330: int DAAppSetQuadraticObjective(TAO_APPLICATION daapplication, PetscTruth pflag){
331:   int info;
332:   TaoMultiGridApplication *mgdaapp;

335:   info=TaoAppGetMultiGridApplication(daapplication,&mgdaapp); CHKERRQ(info);
336:   //  mgdaapp->IsQuadratic=pflag;
337:   if (pflag==PETSC_TRUE){
338:     info = PetscLogInfo((daapplication,"User labels the objective function as quadratic.\n")); CHKERRQ(info);
339:     info = PetscLogInfo((daapplication,"Set KSP MatStructure to SAME_PRECONDITIONER.\n")); CHKERRQ(info);
340:     info=DAAppSetMatStructure(daapplication, SAME_PRECONDITIONER);CHKERRQ(info);
341:   } else {
342:     info = PetscLogInfo((daapplication,"User labels the objective function as NOT quadratic.\n")); CHKERRQ(info);
343:     info = PetscLogInfo((daapplication,"Set KSP MatStructure to SAME_NONZERO_PATTERN.\n")); CHKERRQ(info);
344:     info=DAAppSetMatStructure(daapplication, SAME_NONZERO_PATTERN);CHKERRQ(info);
345:   }
346:   return(0);
347: }


352: /* @
353:   DAAppSetMultiGridOptions - Sets various multigrid options to be used in this application
354:   and the TAO solver.

356:    Collective on TAO_APPLICATION

358:    Input Parameters:
359: .  daapplication - the DA Application object

361:    Level: beginner

363: .keywords:  options

365: .seealso: TaoDAAppSetOptions();

367: @ */
368: int DAAppSetMultiGridOptions(TAO_APPLICATION daapplication){
369:   int info,nlevel;

371:   PetscTruth flg1=PETSC_FALSE;

374:   info = PetscLogInfo((daapplication,"TaoDAAppSetMultiGridOptions(): Reading command line for options\n")); CHKERRQ(info);

376:   flg1=PETSC_FALSE;
377:   info = PetscOptionsInt("-tao_da_multigrid","use multigrid linear solver","DAAppUseMultigrid",
378:                          PETSC_FALSE,&nlevel,&flg1);CHKERRQ(info);
379:   if (flg1) {
380:     info=DAAppUseMultigrid(daapplication,nlevel);CHKERRQ(info);
381:     info = PetscLogInfo((daapplication,"TaoDAAppSetOptions: Use Multigrid linear solver \n")); CHKERRQ(info);
382:   }
383:   /*
384:   flg1=PETSC_FALSE,flg2=PETSC_FALSE;
385:   info = PetscOptionsTruth("-tao_da_quadratic","Identify the objective function as quadratic",
386:                              "DAAppSetQuadraticObjective",PETSC_FALSE,&flg2,&flg1);CHKERRQ(info);
387:   if (flg1) {
388:     info = DAAppSetQuadraticObjective(daapplication, flg2);CHKERRQ(info); 
389:   }
390:   */

392:   return(0);
393: }