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: }