Actual source code: daapp_mgrid.c

  1: #include "tao_general.h"
  2: #include "taodaapplication.h"     /* taodaapplication.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:   PetscInt 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: {
 52:   TaoVecPetsc *px = dynamic_cast <TaoVecPetsc *> (xx);
 53:   TaoMatPetsc *pH = dynamic_cast <TaoMatPetsc *> (HH);

 55:   int     info;
 56:   PetscInt i,clevel,cclevel;
 57:   Vec X = px->GetVec();
 58:   Mat H, HPre;
 59:   MatStructure pflag;
 60:   DA_APPLICATION daapp;

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

 65:   PetscStackPush("TAO User DA Hessian");
 66:   info=pH->GetMatrix(&H,&HPre,&pflag);CHKERRQ(info);
 67:   info=DAAppGetCurrentLevel(this->papp,&clevel); CHKERRQ(info);
 68:   cclevel=PetscMin(clevel,this->coarselevel);
 69:   cclevel=PetscMax(cclevel,0);

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

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

 92:   return(0);
 93: }

 95: int DAAppSetMultigridKSP(GridCtx *, PetscInt, KSP);
 96: int DAAppSetupMultigridMonitor(TAO_APPLICATION, DA, PetscInt, void *);
 97: int TaoAppGetMultiGridApplication(TAO_APPLICATION, TaoMultiGridApplication *);
 98: static char TaoPetscPCMGDAAppXXX[] = "TaoPetscPCMGDAApp";
 99: static char TaoPetscAppXXX[] = "TaoPetscApp";


104: static int DAAppDestroyTaoAppXXX(void*ctx){
105:   TaoMultiGridApplication *mctx=(TaoMultiGridApplication*)ctx;
107:   if (mctx){ delete mctx;}
108:   return(0);
109: }

113: int TaoAppGetMultiGridApplication(TAO_APPLICATION daapplication, TaoMultiGridApplication **mgdaapp){
114:   int info;
115:   PetscInt clevel,ii;
116:   Mat H,HH;
117:   MatStructure         pflag=SAME_NONZERO_PATTERN;
118:   DA_APPLICATION       daapp;
119:   TaoMultiGridApplication*   daapp2;
120:   TaoPetscApplication  *pppappp;

124:   info=TaoAppQueryForObject(daapplication,TaoPetscPCMGDAAppXXX,(void**)&daapp2); CHKERRQ(info);
125:   if (daapp2==0){
126:     daapp2=new TaoMultiGridApplication(PETSC_COMM_WORLD);
127:     daapp2->papp=daapplication;
128:     info=TaoAppDAApp(daapplication,&daapp); CHKERRQ(info);

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

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

156:   Collective on TAO_APPLICATION
157:   
158:   Input Parameters:
159: +  daapplication - the DA Application object
160: -  coarselevel - the coarsest grid to be used in the multigrid preconditioner. (Grid 0 is the coarsest grid.

162:    Level: intermediate

164:    Options Database Key:
165: +  -tao_da_multigrid - use multigrid linear solver
166: -  -ksp_view - view the linear solver

168: Note:
169:   This function should be called after DAAppSetHessianRoutine();

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

174: Note:
175:   Multigrid functionality is still under developement for good performance.

177: .seealso: TaoAppGetKSP(), DAAppSetupMultigrid()

179:    Options Database Key:
180: .  -tao_da_multigrid

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

184: @*/
185: int DAAppUseMultigrid(TAO_APPLICATION daapplication, PetscInt coarselevel){
186:   int info;
188:   UsePCMGGG=coarselevel;
189:   info=DAAppSetBeforeMonitor(daapplication,DAAppSetupMultigridMonitor,0); 
190:   CHKERRQ(info);
191:   return(0);
192: }


197: int DAAppSetupMultigridMonitor(TAO_APPLICATION daapplication, DA da, PetscInt level, void *ctx){
198:   int info;
200:   info = DAAppSetupMultigrid(daapplication,UsePCMGGG); CHKERRQ(info);
201:   return(0);
202: }


207: /*@
208:    DAAppSetupMultigrid - Sets up the multigrid preconditioner.

210:    Collective on TAO_APPLICATION

212:    Input Parameters:
213: +  daapplication - the TAO_APPLICATION context
214: -  coarselevel - the coarsest grid that should be used in the multigrid (>=0)

216:    Note:
217:    Usually the coarselevel is set to 0;

219:    Level: intermediate

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

248:     info=TaoAppGetTaoSolver(daapplication,&tao); CHKERRQ(info);
249:     info=TaoSetupApplicationSolver(daapplication,tao); CHKERRQ(info);
250:     info=TaoAppGetKSP(daapplication,&ksp); CHKERRQ(info);
251:     if (ksp){
253:       info = DAAppSetMultigridKSP(daapp->grid+coarselevel,nn,ksp);CHKERRQ(info);
254:     }
255:     daapp->grid[level].mgrid=1;
256:   }
257:   return(0);
258: }


261:   
264: int DAAppSetMultigridKSP(GridCtx* dagrid, PetscInt nlevels, KSP ksp){
265:   int info;
266:   PetscInt i;
267:   PC pc,pc2;
268:   KSP sksp;
269:   const KSPType ksptype;
270:   const PCType pctype;
271:   PCMGType form;

274:   if (ksp==0 || nlevels<=1){
275:     return(0);
276:   }

279:   info = PetscInfo(ksp,"Set multigrid precondition in optimization solver.\n"); CHKERRQ(info);
280:   
281:   info = KSPGetType(ksp,&ksptype); CHKERRQ(info);
282:   
283:   info = KSPGetPC(ksp,&pc); CHKERRQ(info);
284:   info = PCGetType(pc,&pctype); CHKERRQ(info);
285:   info = PCSetType(pc,PCMG); CHKERRQ(info);
286:   
287:   info = PCMGSetLevels(pc,nlevels,PETSC_NULL); CHKERRQ(info);

289:   for (i=0; i<nlevels; i++) {
290:     
291:     info = PCMGGetSmoother(pc,i,&sksp); CHKERRQ(info);
292:     info = KSPSetType(sksp,ksptype); CHKERRQ(info);

294:     info = KSPGetPC(sksp,&pc2); CHKERRQ(info);
295:     info = PCSetType(pc2,PCJACOBI); CHKERRQ(info);

297:     form=PC_MG_MULTIPLICATIVE;
298:     info=PCMGSetType(pc,form);CHKERRQ(info);

300:     info = KSPSetOperators(sksp,dagrid[i].H,dagrid[i].H,SAME_NONZERO_PATTERN); CHKERRQ(info);
301:     
302:     if (i<nlevels-1){
303:       info = PCMGSetX(pc,i,dagrid[i].X); CHKERRQ(info);
304:       info = PCMGSetRhs(pc,i,dagrid[i].RHS); CHKERRQ(info);      
305:     }
306:     if (i>0){
307:       info = PCMGSetR(pc,i,dagrid[i].R); CHKERRQ(info);
308:     }
309:     info = PCMGSetResidual(pc,i,PCMGDefaultResidual,dagrid[i].H); CHKERRQ(info);
310:   }
311:   for (i=1; i<nlevels; i++) {
312:     info = PCMGSetInterpolation(pc,i,dagrid[i].Interpolate); CHKERRQ(info);
313:     info = PCMGSetRestriction(pc,i,dagrid[i].Interpolate); CHKERRQ(info);
314:   }
315:   
316:   return(0);
317: }

321: /* @
322:   DAAppSetQuadraticObjective - identify the objective function as quadratic or not.

324:   Collective on TAO_APPLICATION

326:   Input Parameters:
327: +  daapplication - the DA Application object
328: -  flag - indicates whether the objective is quadratic or not

330:    Level: intermediate

332:    Options Database Key:
333: .  -tao_da_quadratic - use multigrid linear solver

335:    Note:
336:    The default value is PETSC_FALSE.

338:    Note:
339:    If quadratic, consider setting the flag used in KSP to SAME_PRECONDITIONER

341:    Note:
342:    If quadratic, this routine reduces the number of Hessian evaluations done when using
343:    the multigrid preconditioner.

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


348: .keywords: DA, Objective Function

350: @ */
351: int DAAppSetQuadraticObjective(TAO_APPLICATION daapplication, PetscTruth pflag){
352:   int info;
353:   TaoMultiGridApplication *mgdaapp;

356:   info=TaoAppGetMultiGridApplication(daapplication,&mgdaapp); CHKERRQ(info);
357:   //  mgdaapp->IsQuadratic=pflag;
358:   if (pflag==PETSC_TRUE){
359:     info = PetscInfo(daapplication,"User labels the objective function as quadratic.\n"); CHKERRQ(info);
360:     info = PetscInfo(daapplication,"Set KSP MatStructure to SAME_PRECONDITIONER.\n"); CHKERRQ(info);
361:     info=DAAppSetMatStructure(daapplication, SAME_PRECONDITIONER);CHKERRQ(info);
362:   } else {
363:     info = PetscInfo(daapplication,"User labels the objective function as NOT quadratic.\n"); CHKERRQ(info);
364:     info = PetscInfo(daapplication,"Set KSP MatStructure to SAME_NONZERO_PATTERN.\n"); CHKERRQ(info);
365:     info=DAAppSetMatStructure(daapplication, SAME_NONZERO_PATTERN);CHKERRQ(info);
366:   }
367:   return(0);
368: }


373: /* @
374:   DAAppSetMultiGridOptions - Sets various multigrid options to be used in this application
375:   and the TAO solver.

377:    Collective on TAO_APPLICATION

379:    Input Parameters:
380: .  daapplication - the DA Application object

382:    Level: beginner

384: .keywords:  options

386: .seealso: TaoDAAppSetOptions();

388: @ */
389: int DAAppSetMultiGridOptions(TAO_APPLICATION daapplication){
390:   int info;
391:   PetscInt nlevel;

393:   PetscTruth flg1=PETSC_FALSE;

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

398:   flg1=PETSC_FALSE;
399:   info = PetscOptionsInt("-tao_da_multigrid","use multigrid linear solver","DAAppUseMultigrid",
400:                          PETSC_FALSE,&nlevel,&flg1);CHKERRQ(info);
401:   if (flg1) {
402:     info=DAAppUseMultigrid(daapplication,nlevel);CHKERRQ(info);
403:     info = PetscInfo(daapplication,"TaoDAAppSetOptions: Use Multigrid linear solver \n"); CHKERRQ(info);
404:   }
405:   /*
406:   flg1=PETSC_FALSE,flg2=PETSC_FALSE;
407:   info = PetscOptionsTruth("-tao_da_quadratic","Identify the objective function as quadratic",
408:                              "DAAppSetQuadraticObjective",PETSC_FALSE,&flg2,&flg1);CHKERRQ(info);
409:   if (flg1) {
410:     info = DAAppSetQuadraticObjective(daapplication, flg2);CHKERRQ(info); 
411:   }
412:   */

414:   return(0);
415: }