Actual source code: mg.c

      1: #define PETSCKSP_DLL

      3: /*
      4:     Defines the multigrid preconditioner interface.
      5: */
     6:  #include src/ksp/pc/impls/mg/mgimpl.h

     11: PetscErrorCode PCMGMCycle_Private(PC_MG **mglevels,PetscTruth *converged)
     12: {
     13:   PC_MG          *mg = *mglevels,*mgc;
     15:   PetscInt       cycles = mg->cycles;
     16:   PetscScalar    zero = 0.0;

     19:   if (converged) *converged = PETSC_FALSE;

     21:   if (mg->eventsolve) {PetscLogEventBegin(mg->eventsolve,0,0,0,0);}
     22:   KSPSolve(mg->smoothd,mg->b,mg->x);
     23:   if (mg->eventsolve) {PetscLogEventEnd(mg->eventsolve,0,0,0,0);}
     24:   if (mg->level) {  /* not the coarsest grid */
     25:     (*mg->residual)(mg->A,mg->b,mg->x,mg->r);

     27:     /* if on finest level and have convergence criteria set */
     28:     if (mg->level == mg->levels-1 && mg->ttol) {
     29:       PetscReal rnorm;
     30:       VecNorm(mg->r,NORM_2,&rnorm);
     31:       if (rnorm <= mg->ttol) {
     32:         *converged = PETSC_TRUE;
     33:         if (rnorm < mg->abstol) {
     34:           PetscLogInfo((0,"PCMGMCycle_Private:Linear solver has converged. Residual norm %g is less than absolute tolerance %g\n",rnorm,mg->abstol));
     35:         } else {
     36:           PetscLogInfo((0,"PCMGMCycle_Private:Linear solver has converged. Residual norm %g is less than relative tolerance times initial residual norm %g\n",rnorm,mg->ttol));
     37:         }
     38:         return(0);
     39:       }
     40:     }

     42:     mgc = *(mglevels - 1);
     43:     MatRestrict(mg->restrct,mg->r,mgc->b);
     44:     VecSet(mgc->x,zero);
     45:     while (cycles--) {
     46:       PCMGMCycle_Private(mglevels-1,converged);
     47:     }
     48:     MatInterpolateAdd(mg->interpolate,mgc->x,mg->x,mg->x);
     49:     if (mg->eventsolve) {PetscLogEventBegin(mg->eventsolve,0,0,0,0);}
     50:     KSPSolve(mg->smoothu,mg->b,mg->x);
     51:     if (mg->eventsolve) {PetscLogEventEnd(mg->eventsolve,0,0,0,0);}
     52:   }
     53:   return(0);
     54: }

     56: /*
     57:        PCMGCreate_Private - Creates a PC_MG structure for use with the
     58:                multigrid code. Level 0 is the coarsest. (But the
     59:                finest level is stored first in the array).

     61: */
     64: static PetscErrorCode PCMGCreate_Private(MPI_Comm comm,PetscInt levels,PC pc,MPI_Comm *comms,PC_MG ***result)
     65: {
     66:   PC_MG          **mg;
     68:   PetscInt       i;
     69:   PetscMPIInt    size;
     70:   const char     *prefix;
     71:   PC             ipc;

     74:   PetscMalloc(levels*sizeof(PC_MG*),&mg);
     75:   PetscLogObjectMemory(pc,levels*(sizeof(PC_MG*)+sizeof(PC_MG)));

     77:   PCGetOptionsPrefix(pc,&prefix);

     79:   for (i=0; i<levels; i++) {
     80:     PetscNew(PC_MG,&mg[i]);
     81:     mg[i]->level           = i;
     82:     mg[i]->levels          = levels;
     83:     mg[i]->cycles          = 1;
     84:     mg[i]->galerkin        = PETSC_FALSE;
     85:     mg[i]->galerkinused    = PETSC_FALSE;
     86:     mg[i]->default_smoothu = 1;
     87:     mg[i]->default_smoothd = 1;

     89:     if (comms) comm = comms[i];
     90:     KSPCreate(comm,&mg[i]->smoothd);
     91:     KSPSetTolerances(mg[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT, mg[i]->default_smoothd);
     92:     KSPSetOptionsPrefix(mg[i]->smoothd,prefix);

     94:     /* do special stuff for coarse grid */
     95:     if (!i && levels > 1) {
     96:       KSPAppendOptionsPrefix(mg[0]->smoothd,"mg_coarse_");

     98:       /* coarse solve is (redundant) LU by default */
     99:       KSPSetType(mg[0]->smoothd,KSPPREONLY);
    100:       KSPGetPC(mg[0]->smoothd,&ipc);
    101:       MPI_Comm_size(comm,&size);
    102:       if (size > 1) {
    103:         PCSetType(ipc,PCREDUNDANT);
    104:         PCRedundantGetPC(ipc,&ipc);
    105:       }
    106:       PCSetType(ipc,PCLU);

    108:     } else {
    109:       char tprefix[128];
    110:       sprintf(tprefix,"mg_levels_%d_",(int)i);
    111:       KSPAppendOptionsPrefix(mg[i]->smoothd,tprefix);
    112:     }
    113:     PetscLogObjectParent(pc,mg[i]->smoothd);
    114:     mg[i]->smoothu         = mg[i]->smoothd;
    115:     mg[i]->rtol = 0.0;
    116:     mg[i]->abstol = 0.0;
    117:     mg[i]->dtol = 0.0;
    118:     mg[i]->ttol = 0.0;
    119:     mg[i]->eventsetup = 0;
    120:     mg[i]->eventsolve = 0;
    121:   }
    122:   *result = mg;
    123:   return(0);
    124: }

    128: static PetscErrorCode PCDestroy_MG(PC pc)
    129: {
    130:   PC_MG          **mg = (PC_MG**)pc->data;
    132:   PetscInt       i,n = mg[0]->levels;

    135:   if (mg[0]->galerkinused) {
    136:     Mat B;
    137:     for (i=0; i<n-1; i++) {
    138:       KSPGetOperators(mg[i]->smoothd,0,&B,0);
    139:       MatDestroy(B);
    140:     }
    141:   }

    143:   for (i=0; i<n-1; i++) {
    144:     if (mg[i+1]->r) {VecDestroy(mg[i+1]->r);}
    145:     if (mg[i]->b) {VecDestroy(mg[i]->b);}
    146:     if (mg[i]->x) {VecDestroy(mg[i]->x);}
    147:     if (mg[i+1]->restrct) {MatDestroy(mg[i+1]->restrct);}
    148:     if (mg[i+1]->interpolate) {MatDestroy(mg[i+1]->interpolate);}
    149:   }

    151:   for (i=0; i<n; i++) {
    152:     if (mg[i]->smoothd != mg[i]->smoothu) {
    153:       KSPDestroy(mg[i]->smoothd);
    154:     }
    155:     KSPDestroy(mg[i]->smoothu);
    156:     PetscFree(mg[i]);
    157:   }
    158:   PetscFree(mg);
    159:   return(0);
    160: }

    164: EXTERN PetscErrorCode PCMGACycle_Private(PC_MG**);
    165: EXTERN PetscErrorCode PCMGFCycle_Private(PC_MG**);
    166: EXTERN PetscErrorCode PCMGKCycle_Private(PC_MG**);

    168: /*
    169:    PCApply_MG - Runs either an additive, multiplicative, Kaskadic
    170:              or full cycle of multigrid.

    172:   Note:
    173:   A simple wrapper which calls PCMGMCycle(),PCMGACycle(), or PCMGFCycle().
    174: */
    177: static PetscErrorCode PCApply_MG(PC pc,Vec b,Vec x)
    178: {
    179:   PC_MG          **mg = (PC_MG**)pc->data;
    180:   PetscScalar    zero = 0.0;
    182:   PetscInt       levels = mg[0]->levels;

    185:   mg[levels-1]->b = b;
    186:   mg[levels-1]->x = x;
    187:   if (!mg[levels-1]->r && mg[0]->am == PC_MG_ADDITIVE) {
    188:     Vec tvec;
    189:     VecDuplicate(mg[levels-1]->b,&tvec);
    190:     PCMGSetR(pc,levels-1,tvec);
    191:     VecDestroy(tvec);
    192:   }
    193:   if (mg[0]->am == PC_MG_MULTIPLICATIVE) {
    194:     VecSet(x,zero);
    195:     PCMGMCycle_Private(mg+levels-1,PETSC_NULL);
    196:   }
    197:   else if (mg[0]->am == PC_MG_ADDITIVE) {
    198:     PCMGACycle_Private(mg);
    199:   }
    200:   else if (mg[0]->am == PC_MG_KASKADE) {
    201:     PCMGKCycle_Private(mg);
    202:   }
    203:   else {
    204:     PCMGFCycle_Private(mg);
    205:   }
    206:   return(0);
    207: }

    211: static PetscErrorCode PCApplyRichardson_MG(PC pc,Vec b,Vec x,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its)
    212: {
    213:   PC_MG          **mg = (PC_MG**)pc->data;
    215:   PetscInt       levels = mg[0]->levels;
    216:   PetscTruth     converged = PETSC_FALSE;

    219:   mg[levels-1]->b    = b;
    220:   mg[levels-1]->x    = x;

    222:   mg[levels-1]->rtol = rtol;
    223:   mg[levels-1]->abstol = abstol;
    224:   mg[levels-1]->dtol = dtol;
    225:   if (rtol) {
    226:     /* compute initial residual norm for relative convergence test */
    227:     PetscReal rnorm;
    228:     (*mg[levels-1]->residual)(mg[levels-1]->A,b,x,w);
    229:     VecNorm(w,NORM_2,&rnorm);
    230:     mg[levels-1]->ttol = PetscMax(rtol*rnorm,abstol);
    231:   } else if (abstol) {
    232:     mg[levels-1]->ttol = abstol;
    233:   } else {
    234:     mg[levels-1]->ttol = 0.0;
    235:   }

    237:   while (its-- && !converged) {
    238:     PCMGMCycle_Private(mg+levels-1,&converged);
    239:   }
    240:   return(0);
    241: }

    245: PetscErrorCode PCSetFromOptions_MG(PC pc)
    246: {
    248:   PetscInt       m,levels = 1;
    249:   PetscTruth     flg;
    250:   PC_MG          **mg = (PC_MG**)pc->data;
    251:   PCMGType       mgtype = mg[0]->am;;

    255:   PetscOptionsHead("Multigrid options");
    256:     if (!pc->data) {
    257:       PetscOptionsInt("-pc_mg_levels","Number of Levels","PCMGSetLevels",levels,&levels,&flg);
    258:       PCMGSetLevels(pc,levels,PETSC_NULL);
    259:     }
    260:     PetscOptionsInt("-pc_mg_cycles","1 for V cycle, 2 for W-cycle","PCMGSetCycles",1,&m,&flg);
    261:     if (flg) {
    262:       PCMGSetCycles(pc,m);
    263:     }
    264:     PetscOptionsName("-pc_mg_galerkin","Use Galerkin process to compute coarser operators","PCMGSetGalerkin",&flg);
    265:     if (flg) {
    266:       PCMGSetGalerkin(pc);
    267:     }
    268:     PetscOptionsInt("-pc_mg_smoothup","Number of post-smoothing steps","PCMGSetNumberSmoothUp",1,&m,&flg);
    269:     if (flg) {
    270:       PCMGSetNumberSmoothUp(pc,m);
    271:     }
    272:     PetscOptionsInt("-pc_mg_smoothdown","Number of pre-smoothing steps","PCMGSetNumberSmoothDown",1,&m,&flg);
    273:     if (flg) {
    274:       PCMGSetNumberSmoothDown(pc,m);
    275:     }
    276:     PetscOptionsEnum("-pc_mg_type","Multigrid type","PCMGSetType",PCMGTypes,(PetscEnum)mgtype,(PetscEnum*)&mgtype,&flg);
    277:     if (flg) {PCMGSetType(pc,mgtype);}
    278:     PetscOptionsName("-pc_mg_log","Log times for each multigrid level","None",&flg);
    279:     if (flg) {
    280:       PetscInt i;
    281:       char     eventname[128];
    282:       if (!mg) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
    283:       levels = mg[0]->levels;
    284:       for (i=0; i<levels; i++) {
    285:         sprintf(eventname,"MSetup Level %d",(int)i);
    286:         PetscLogEventRegister(&mg[i]->eventsetup,eventname,pc->cookie);
    287:         sprintf(eventname,"MGSolve Level %d to 0",(int)i);
    288:         PetscLogEventRegister(&mg[i]->eventsolve,eventname,pc->cookie);
    289:       }
    290:     }
    291:   PetscOptionsTail();
    292:   return(0);
    293: }

    295: const char *PCMGTypes[] = {"MULTIPLICATIVE","ADDITIVE","FULL","KASKADE","PCMGType","PC_MG",0};

    299: static PetscErrorCode PCView_MG(PC pc,PetscViewer viewer)
    300: {
    301:   PC_MG          **mg = (PC_MG**)pc->data;
    303:   PetscInt       levels = mg[0]->levels,i;
    304:   PetscTruth     iascii;

    307:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
    308:   if (iascii) {
    309:     PetscViewerASCIIPrintf(viewer,"  MG: type is %s, levels=%D cycles=%D, pre-smooths=%D, post-smooths=%D\n",
    310:                       PCMGTypes[mg[0]->am],levels,mg[0]->cycles,mg[0]->default_smoothd,mg[0]->default_smoothu);
    311:     if (mg[0]->galerkin) {
    312:       PetscViewerASCIIPrintf(viewer,"    Using Galerkin computed coarse grid matrices\n");
    313:     }
    314:     for (i=0; i<levels; i++) {
    315:       if (!i) {
    316:         PetscViewerASCIIPrintf(viewer,"Coarse gride solver -- level %D -------------------------------\n",i);
    317:       } else {
    318:         PetscViewerASCIIPrintf(viewer,"Down solver (pre-smoother) on level %D -------------------------------\n",i);
    319:       }
    320:       PetscViewerASCIIPushTab(viewer);
    321:       KSPView(mg[i]->smoothd,viewer);
    322:       PetscViewerASCIIPopTab(viewer);
    323:       if (i && mg[i]->smoothd == mg[i]->smoothu) {
    324:         PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) same as down solver (pre-smoother)\n");
    325:       } else if (i){
    326:         PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) on level %D -------------------------------\n",i);
    327:         PetscViewerASCIIPushTab(viewer);
    328:         KSPView(mg[i]->smoothu,viewer);
    329:         PetscViewerASCIIPopTab(viewer);
    330:       }
    331:     }
    332:   } else {
    333:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCMG",((PetscObject)viewer)->type_name);
    334:   }
    335:   return(0);
    336: }

    338: /*
    339:     Calls setup for the KSP on each level
    340: */
    343: static PetscErrorCode PCSetUp_MG(PC pc)
    344: {
    345:   PC_MG          **mg = (PC_MG**)pc->data;
    347:   PetscInt       i,n = mg[0]->levels;
    348:   PC             cpc;
    349:   PetscTruth     preonly,lu,redundant,cholesky,monitor = PETSC_FALSE,dump;
    350:   PetscViewer    ascii;
    351:   MPI_Comm       comm;
    352:   Mat            dA,dB;
    353:   MatStructure   uflag;
    354:   Vec            tvec;

    357:   if (!pc->setupcalled) {
    358:     PetscOptionsHasName(0,"-pc_mg_monitor",&monitor);
    359:
    360:     for (i=0; i<n; i++) {
    361:       if (monitor) {
    362:         PetscObjectGetComm((PetscObject)mg[i]->smoothd,&comm);
    363:         PetscViewerASCIIOpen(comm,"stdout",&ascii);
    364:         PetscViewerASCIISetTab(ascii,n-i);
    365:         KSPSetMonitor(mg[i]->smoothd,KSPDefaultMonitor,ascii,(PetscErrorCode(*)(void*))PetscViewerDestroy);
    366:       }
    367:       KSPSetFromOptions(mg[i]->smoothd);
    368:     }
    369:     for (i=1; i<n; i++) {
    370:       if (mg[i]->smoothu && mg[i]->smoothu != mg[i]->smoothd) {
    371:         if (monitor) {
    372:           PetscObjectGetComm((PetscObject)mg[i]->smoothu,&comm);
    373:           PetscViewerASCIIOpen(comm,"stdout",&ascii);
    374:           PetscViewerASCIISetTab(ascii,n-i);
    375:           KSPSetMonitor(mg[i]->smoothu,KSPDefaultMonitor,ascii,(PetscErrorCode(*)(void*))PetscViewerDestroy);
    376:         }
    377:         KSPSetFromOptions(mg[i]->smoothu);
    378:       }
    379:     }
    380:     for (i=1; i<n; i++) {
    381:       if (mg[i]->restrct && !mg[i]->interpolate) {
    382:         PCMGSetInterpolate(pc,i,mg[i]->restrct);
    383:       }
    384:       if (!mg[i]->restrct && mg[i]->interpolate) {
    385:         PCMGSetRestriction(pc,i,mg[i]->interpolate);
    386:       }
    387: #if defined(PETSC_USE_DEBUG)
    388:       if (!mg[i]->restrct || !mg[i]->interpolate) {
    389:         SETERRQ1(PETSC_ERR_ARG_WRONGSTATE,"Need to set restriction or interpolation on level %d",(int)i);
    390:       }
    391: #endif
    392:     }
    393:     for (i=0; i<n-1; i++) {
    394:       if (!mg[i]->r && i) {
    395:         VecDuplicate(mg[i]->b,&tvec);
    396:         PCMGSetR(pc,i,tvec);
    397:         VecDestroy(tvec);
    398:       }
    399:       if (!mg[i]->x) {
    400:         VecDuplicate(mg[i]->b,&tvec);
    401:         PCMGSetX(pc,i,tvec);
    402:         VecDestroy(tvec);
    403:       }
    404:     }
    405:   }

    407:   /* If user did not provide fine grid operators, use those from PC */
    408:   /* BUG BUG BUG This will work ONLY the first time called: hence if the user changes
    409:      the PC matrices between solves PCMG will continue to use first set provided */
    410:   KSPGetOperators(mg[n-1]->smoothd,&dA,&dB,&uflag);
    411:   if (!dA  && !dB) {
    412:     PetscLogInfo((pc,"PCSetUp_MG: Using outer operators to define finest grid operator \n  because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n"));
    413:     KSPSetOperators(mg[n-1]->smoothd,pc->mat,pc->pmat,uflag);
    414:   }

    416:   if (mg[0]->galerkin) {
    417:     Mat B;
    418:     mg[0]->galerkinused = PETSC_TRUE;
    419:     /* currently only handle case where mat and pmat are the same on coarser levels */
    420:     KSPGetOperators(mg[n-1]->smoothd,&dA,&dB,&uflag);
    421:     if (!pc->setupcalled) {
    422:       for (i=n-2; i>-1; i--) {
    423:         MatPtAP(dB,mg[i+1]->interpolate,MAT_INITIAL_MATRIX,1.0,&B);
    424:         KSPSetOperators(mg[i]->smoothd,B,B,uflag);
    425:         dB   = B;
    426:       }
    427:     } else {
    428:       for (i=n-2; i>-1; i--) {
    429:         KSPGetOperators(mg[i]->smoothd,0,&B,0);
    430:         MatPtAP(dB,mg[i]->interpolate,MAT_REUSE_MATRIX,1.0,&B);
    431:         KSPSetOperators(mg[i]->smoothd,B,B,uflag);
    432:         dB   = B;
    433:       }
    434:     }
    435:   }

    437:   for (i=1; i<n; i++) {
    438:     if (mg[i]->smoothu == mg[i]->smoothd) {
    439:       /* if doing only down then initial guess is zero */
    440:       KSPSetInitialGuessNonzero(mg[i]->smoothd,PETSC_TRUE);
    441:     }
    442:     if (mg[i]->eventsetup) {PetscLogEventBegin(mg[i]->eventsetup,0,0,0,0);}
    443:     KSPSetUp(mg[i]->smoothd);
    444:     if (mg[i]->eventsetup) {PetscLogEventEnd(mg[i]->eventsetup,0,0,0,0);}
    445:   }
    446:   for (i=1; i<n; i++) {
    447:     if (mg[i]->smoothu && mg[i]->smoothu != mg[i]->smoothd) {
    448:       PC           uppc,downpc;
    449:       Mat          downmat,downpmat,upmat,uppmat;
    450:       MatStructure matflag;

    452:       /* check if operators have been set for up, if not use down operators to set them */
    453:       KSPGetPC(mg[i]->smoothu,&uppc);
    454:       PCGetOperators(uppc,&upmat,&uppmat,PETSC_NULL);
    455:       if (!upmat) {
    456:         KSPGetPC(mg[i]->smoothd,&downpc);
    457:         PCGetOperators(downpc,&downmat,&downpmat,&matflag);
    458:         KSPSetOperators(mg[i]->smoothu,downmat,downpmat,matflag);
    459:       }

    461:       KSPSetInitialGuessNonzero(mg[i]->smoothu,PETSC_TRUE);
    462:       if (mg[i]->eventsetup) {PetscLogEventBegin(mg[i]->eventsetup,0,0,0,0);}
    463:       KSPSetUp(mg[i]->smoothu);
    464:       if (mg[i]->eventsetup) {PetscLogEventEnd(mg[i]->eventsetup,0,0,0,0);}
    465:     }
    466:   }

    468:   /*
    469:       If coarse solver is not direct method then DO NOT USE preonly
    470:   */
    471:   PetscTypeCompare((PetscObject)mg[0]->smoothd,KSPPREONLY,&preonly);
    472:   if (preonly) {
    473:     KSPGetPC(mg[0]->smoothd,&cpc);
    474:     PetscTypeCompare((PetscObject)cpc,PCLU,&lu);
    475:     PetscTypeCompare((PetscObject)cpc,PCREDUNDANT,&redundant);
    476:     PetscTypeCompare((PetscObject)cpc,PCCHOLESKY,&cholesky);
    477:     if (!lu && !redundant && !cholesky) {
    478:       KSPSetType(mg[0]->smoothd,KSPGMRES);
    479:     }
    480:   }

    482:   if (!pc->setupcalled) {
    483:     if (monitor) {
    484:       PetscObjectGetComm((PetscObject)mg[0]->smoothd,&comm);
    485:       PetscViewerASCIIOpen(comm,"stdout",&ascii);
    486:       PetscViewerASCIISetTab(ascii,n);
    487:       KSPSetMonitor(mg[0]->smoothd,KSPDefaultMonitor,ascii,(PetscErrorCode(*)(void*))PetscViewerDestroy);
    488:     }
    489:     KSPSetFromOptions(mg[0]->smoothd);
    490:   }

    492:   if (mg[0]->eventsetup) {PetscLogEventBegin(mg[0]->eventsetup,0,0,0,0);}
    493:   KSPSetUp(mg[0]->smoothd);
    494:   if (mg[0]->eventsetup) {PetscLogEventEnd(mg[0]->eventsetup,0,0,0,0);}

    496: #if defined(PETSC_USE_SOCKET_VIEWER)
    497:   /*
    498:      Dump the interpolation/restriction matrices to matlab plus the
    499:    Jacobian/stiffness on each level. This allows Matlab users to
    500:    easily check if the Galerkin condition A_c = R A_f R^T is satisfied */
    501:   PetscOptionsHasName(pc->prefix,"-pc_mg_dump_matlab",&dump);
    502:   if (dump) {
    503:     for (i=1; i<n; i++) {
    504:       MatView(mg[i]->restrct,PETSC_VIEWER_SOCKET_(pc->comm));
    505:     }
    506:     for (i=0; i<n; i++) {
    507:       KSPGetPC(mg[i]->smoothd,&pc);
    508:       MatView(pc->mat,PETSC_VIEWER_SOCKET_(pc->comm));
    509:     }
    510:   }
    511: #endif

    513:   PetscOptionsHasName(pc->prefix,"-pc_mg_dump_binary",&dump);
    514:   if (dump) {
    515:     for (i=1; i<n; i++) {
    516:       MatView(mg[i]->restrct,PETSC_VIEWER_BINARY_(pc->comm));
    517:     }
    518:     for (i=0; i<n; i++) {
    519:       KSPGetPC(mg[i]->smoothd,&pc);
    520:       MatView(pc->mat,PETSC_VIEWER_BINARY_(pc->comm));
    521:     }
    522:   }
    523:   return(0);
    524: }

    526: /* -------------------------------------------------------------------------------------*/

    530: /*@C
    531:    PCMGSetLevels - Sets the number of levels to use with MG.
    532:    Must be called before any other MG routine.

    534:    Collective on PC

    536:    Input Parameters:
    537: +  pc - the preconditioner context
    538: .  levels - the number of levels
    539: -  comms - optional communicators for each level; this is to allow solving the coarser problems
    540:            on smaller sets of processors. Use PETSC_NULL_OBJECT for default in Fortran

    542:    Level: intermediate

    544:    Notes:
    545:      If the number of levels is one then the multigrid uses the -mg_levels prefix
    546:   for setting the level options rather than the -mg_coarse prefix.

    548: .keywords: MG, set, levels, multigrid

    550: .seealso: PCMGSetType(), PCMGGetLevels()
    551: @*/
    552: PetscErrorCode PETSCKSP_DLLEXPORT PCMGSetLevels(PC pc,PetscInt levels,MPI_Comm *comms)
    553: {
    555:   PC_MG          **mg;

    560:   if (pc->data) {
    561:     SETERRQ(PETSC_ERR_ORDER,"Number levels already set for MG\n\
    562:     make sure that you call PCMGSetLevels() before KSPSetFromOptions()");
    563:   }
    564:   PCMGCreate_Private(pc->comm,levels,pc,comms,&mg);
    565:   mg[0]->am                = PC_MG_MULTIPLICATIVE;
    566:   pc->data                 = (void*)mg;
    567:   pc->ops->applyrichardson = PCApplyRichardson_MG;
    568:   return(0);
    569: }

    573: /*@
    574:    PCMGGetLevels - Gets the number of levels to use with MG.

    576:    Not Collective

    578:    Input Parameter:
    579: .  pc - the preconditioner context

    581:    Output parameter:
    582: .  levels - the number of levels

    584:    Level: advanced

    586: .keywords: MG, get, levels, multigrid

    588: .seealso: PCMGSetLevels()
    589: @*/
    590: PetscErrorCode PETSCKSP_DLLEXPORT PCMGGetLevels(PC pc,PetscInt *levels)
    591: {
    592:   PC_MG  **mg;

    598:   mg      = (PC_MG**)pc->data;
    599:   *levels = mg[0]->levels;
    600:   return(0);
    601: }

    605: /*@
    606:    PCMGSetType - Determines the form of multigrid to use:
    607:    multiplicative, additive, full, or the Kaskade algorithm.

    609:    Collective on PC

    611:    Input Parameters:
    612: +  pc - the preconditioner context
    613: -  form - multigrid form, one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE,
    614:    PC_MG_FULL, PC_MG_KASKADE

    616:    Options Database Key:
    617: .  -pc_mg_type <form> - Sets <form>, one of multiplicative,
    618:    additive, full, kaskade

    620:    Level: advanced

    622: .keywords: MG, set, method, multiplicative, additive, full, Kaskade, multigrid

    624: .seealso: PCMGSetLevels()
    625: @*/
    626: PetscErrorCode PETSCKSP_DLLEXPORT PCMGSetType(PC pc,PCMGType form)
    627: {
    628:   PC_MG **mg;

    632:   mg = (PC_MG**)pc->data;

    634:   if (!mg) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
    635:   mg[0]->am = form;
    636:   if (form == PC_MG_MULTIPLICATIVE) pc->ops->applyrichardson = PCApplyRichardson_MG;
    637:   else pc->ops->applyrichardson = 0;
    638:   return(0);
    639: }

    643: /*@
    644:    PCMGSetCycles - Sets the type cycles to use.  Use PCMGSetCyclesOnLevel() for more
    645:    complicated cycling.

    647:    Collective on PC

    649:    Input Parameters:
    650: +  pc - the multigrid context
    651: -  n - the number of cycles

    653:    Options Database Key:
    654: $  -pc_mg_cycles n - 1 denotes a V-cycle; 2 denotes a W-cycle.

    656:    Level: advanced

    658: .keywords: MG, set, cycles, V-cycle, W-cycle, multigrid

    660: .seealso: PCMGSetCyclesOnLevel()
    661: @*/
    662: PetscErrorCode PETSCKSP_DLLEXPORT PCMGSetCycles(PC pc,PetscInt n)
    663: {
    664:   PC_MG    **mg;
    665:   PetscInt i,levels;

    669:   mg     = (PC_MG**)pc->data;
    670:   if (!mg) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
    671:   levels = mg[0]->levels;

    673:   for (i=0; i<levels; i++) {
    674:     mg[i]->cycles  = n;
    675:   }
    676:   return(0);
    677: }

    681: /*@
    682:    PCMGSetGalerkin - Causes the coarser grid matrices to be computed from the
    683:       finest grid via the Galerkin process: A_i-1 = r_i * A_i * r_i^t

    685:    Collective on PC

    687:    Input Parameters:
    688: +  pc - the multigrid context
    689: -  n - the number of cycles

    691:    Options Database Key:
    692: $  -pc_mg_galerkin

    694:    Level: intermediate

    696: .keywords: MG, set, Galerkin

    698: @*/
    699: PetscErrorCode PETSCKSP_DLLEXPORT PCMGSetGalerkin(PC pc)
    700: {
    701:   PC_MG    **mg;
    702:   PetscInt i,levels;

    706:   mg     = (PC_MG**)pc->data;
    707:   if (!mg) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
    708:   levels = mg[0]->levels;

    710:   for (i=0; i<levels; i++) {
    711:     mg[i]->galerkin = PETSC_TRUE;
    712:   }
    713:   return(0);
    714: }

    718: /*@
    719:    PCMGSetNumberSmoothDown - Sets the number of pre-smoothing steps to
    720:    use on all levels. Use PCMGGetSmootherDown() to set different
    721:    pre-smoothing steps on different levels.

    723:    Collective on PC

    725:    Input Parameters:
    726: +  mg - the multigrid context
    727: -  n - the number of smoothing steps

    729:    Options Database Key:
    730: .  -pc_mg_smoothdown <n> - Sets number of pre-smoothing steps

    732:    Level: advanced

    734: .keywords: MG, smooth, down, pre-smoothing, steps, multigrid

    736: .seealso: PCMGSetNumberSmoothUp()
    737: @*/
    738: PetscErrorCode PETSCKSP_DLLEXPORT PCMGSetNumberSmoothDown(PC pc,PetscInt n)
    739: {
    740:   PC_MG          **mg;
    742:   PetscInt       i,levels;

    746:   mg     = (PC_MG**)pc->data;
    747:   if (!mg) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
    748:   levels = mg[0]->levels;

    750:   for (i=0; i<levels; i++) {
    751:     /* make sure smoother up and down are different */
    752:     PCMGGetSmootherUp(pc,i,PETSC_NULL);
    753:     KSPSetTolerances(mg[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);
    754:     mg[i]->default_smoothd = n;
    755:   }
    756:   return(0);
    757: }

    761: /*@
    762:    PCMGSetNumberSmoothUp - Sets the number of post-smoothing steps to use
    763:    on all levels. Use PCMGGetSmootherUp() to set different numbers of
    764:    post-smoothing steps on different levels.

    766:    Collective on PC

    768:    Input Parameters:
    769: +  mg - the multigrid context
    770: -  n - the number of smoothing steps

    772:    Options Database Key:
    773: .  -pc_mg_smoothup <n> - Sets number of post-smoothing steps

    775:    Level: advanced

    777:    Note: this does not set a value on the coarsest grid, since we assume that
    778:     there is no seperate smooth up on the coarsest grid.

    780: .keywords: MG, smooth, up, post-smoothing, steps, multigrid

    782: .seealso: PCMGSetNumberSmoothDown()
    783: @*/
    784: PetscErrorCode PETSCKSP_DLLEXPORT PCMGSetNumberSmoothUp(PC pc,PetscInt n)
    785: {
    786:   PC_MG          **mg;
    788:   PetscInt       i,levels;

    792:   mg     = (PC_MG**)pc->data;
    793:   if (!mg) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
    794:   levels = mg[0]->levels;

    796:   for (i=1; i<levels; i++) {
    797:     /* make sure smoother up and down are different */
    798:     PCMGGetSmootherUp(pc,i,PETSC_NULL);
    799:     KSPSetTolerances(mg[i]->smoothu,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);
    800:     mg[i]->default_smoothu = n;
    801:   }
    802:   return(0);
    803: }

    805: /* ----------------------------------------------------------------------------------------*/

    807: /*MC
    808:    PCMG - Use geometric multigrid preconditioning. This preconditioner requires you provide additional
    809:     information about the coarser grid matrices and restriction/interpolation operators.

    811:    Options Database Keys:
    812: +  -pc_mg_levels <nlevels> - number of levels including finest
    813: .  -pc_mg_cycles 1 or 2 - for V or W-cycle
    814: .  -pc_mg_smoothup <n> - number of smoothing steps after interpolation
    815: .  -pc_mg_smoothdown <n> - number of smoothing steps before applying restriction operator
    816: .  -pc_mg_type <additive,multiplicative,full,cascade> - multiplicative is the default
    817: .  -pc_mg_log - log information about time spent on each level of the solver
    818: .  -pc_mg_monitor - print information on the multigrid convergence
    819: .  -pc_mg_galerkin - use Galerkin process to compute coarser operators
    820: -  -pc_mg_dump_matlab - dumps the matrices for each level and the restriction/interpolation matrices
    821:                         to the Socket viewer for reading from Matlab.

    823:    Notes:

    825:    Level: intermediate

    827:    Concepts: multigrid

    829: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, PCMGType,
    830:            PCMGSetLevels(), PCMGGetLevels(), PCMGSetType(), PCMGSetCycles(), PCMGSetNumberSmoothDown(),
    831:            PCMGSetNumberSmoothUp(), PCMGGetCoarseSolve(), PCMGSetResidual(), PCMGSetInterpolation(),
    832:            PCMGSetRestriction(), PCMGGetSmoother(), PCMGGetSmootherUp(), PCMGGetSmootherDown(),
    833:            PCMGSetCyclesOnLevel(), PCMGSetRhs(), PCMGSetX(), PCMGSetR()
    834: M*/

    839: PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_MG(PC pc)
    840: {
    842:   pc->ops->apply          = PCApply_MG;
    843:   pc->ops->setup          = PCSetUp_MG;
    844:   pc->ops->destroy        = PCDestroy_MG;
    845:   pc->ops->setfromoptions = PCSetFromOptions_MG;
    846:   pc->ops->view           = PCView_MG;

    848:   pc->data                = (void*)0;
    849:   return(0);
    850: }