Actual source code: sor.c

      1: #define PETSCKSP_DLL

      3: /*
      4:    Defines a  (S)SOR  preconditioner for any Mat implementation
      5: */
     6:  #include src/ksp/pc/pcimpl.h

      8: typedef struct {
      9:   PetscInt   its;        /* inner iterations, number of sweeps */
     10:   PetscInt   lits;       /* local inner iterations, number of sweeps applied by the local matrix mat->A */
     11:   MatSORType sym;        /* forward, reverse, symmetric etc. */
     12:   PetscReal  omega;
     13: } PC_SOR;

     17: static PetscErrorCode PCDestroy_SOR(PC pc)
     18: {
     19:   PC_SOR         *jac = (PC_SOR*)pc->data;

     23:   PetscFree(jac);
     24:   return(0);
     25: }

     29: static PetscErrorCode PCApply_SOR(PC pc,Vec x,Vec y)
     30: {
     31:   PC_SOR         *jac = (PC_SOR*)pc->data;
     33:   PetscInt       flag = jac->sym | SOR_ZERO_INITIAL_GUESS;

     36:   MatRelax(pc->pmat,x,jac->omega,(MatSORType)flag,0.0,jac->its,jac->lits,y);
     37:   return(0);
     38: }

     42: static PetscErrorCode PCApplyRichardson_SOR(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its)
     43: {
     44:   PC_SOR         *jac = (PC_SOR*)pc->data;

     48:   PetscLogInfo((pc,"PCApplyRichardson_SOR: Warning, convergence critera ignored, using %D iterations\n",its));
     49:   its  = its*jac->its;
     50:   MatRelax(pc->pmat,b,jac->omega,(MatSORType)jac->sym,0.0,its,jac->lits,y);
     51:   return(0);
     52: }

     56: PetscErrorCode PCSetFromOptions_SOR(PC pc)
     57: {
     58:   PC_SOR         *jac = (PC_SOR*)pc->data;
     60:   PetscTruth     flg;

     63:   PetscOptionsHead("(S)SOR options");
     64:     PetscOptionsReal("-pc_sor_omega","relaxation factor (0 < omega < 2)","PCSORSetOmega",jac->omega,&jac->omega,0);
     65:     PetscOptionsInt("-pc_sor_its","number of inner SOR iterations","PCSORSetIterations",jac->its,&jac->its,0);
     66:     PetscOptionsInt("-pc_sor_lits","number of local inner SOR iterations","PCSORSetIterations",jac->lits,&jac->lits,0);
     67:     PetscOptionsTruthGroupBegin("-pc_sor_symmetric","SSOR, not SOR","PCSORSetSymmetric",&flg);
     68:     if (flg) {PCSORSetSymmetric(pc,SOR_SYMMETRIC_SWEEP);}
     69:     PetscOptionsTruthGroup("-pc_sor_backward","use backward sweep instead of forward","PCSORSetSymmetric",&flg);
     70:     if (flg) {PCSORSetSymmetric(pc,SOR_BACKWARD_SWEEP);}
     71:     PetscOptionsTruthGroup("-pc_sor_local_symmetric","use SSOR seperately on each processor","PCSORSetSymmetric",&flg);
     72:     if (flg) {PCSORSetSymmetric(pc,SOR_LOCAL_SYMMETRIC_SWEEP);}
     73:     PetscOptionsTruthGroup("-pc_sor_local_backward","use backward sweep locally","PCSORSetSymmetric",&flg);
     74:     if (flg) {PCSORSetSymmetric(pc,SOR_LOCAL_BACKWARD_SWEEP);}
     75:     PetscOptionsTruthGroupEnd("-pc_sor_local_forward","use forward sweep locally","PCSORSetSymmetric",&flg);
     76:     if (flg) {PCSORSetSymmetric(pc,SOR_LOCAL_FORWARD_SWEEP);}
     77:   PetscOptionsTail();
     78:   return(0);
     79: }

     83: PetscErrorCode PCView_SOR(PC pc,PetscViewer viewer)
     84: {
     85:   PC_SOR         *jac = (PC_SOR*)pc->data;
     86:   MatSORType     sym = jac->sym;
     87:   const char     *sortype;
     89:   PetscTruth     iascii;

     92:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
     93:   if (iascii) {
     94:     if (sym & SOR_ZERO_INITIAL_GUESS) {PetscViewerASCIIPrintf(viewer,"  SOR:  zero initial guess\n");}
     95:     if (sym == SOR_APPLY_UPPER)              sortype = "apply_upper";
     96:     else if (sym == SOR_APPLY_LOWER)         sortype = "apply_lower";
     97:     else if (sym & SOR_EISENSTAT)            sortype = "Eisenstat";
     98:     else if ((sym & SOR_SYMMETRIC_SWEEP) == SOR_SYMMETRIC_SWEEP)
     99:                                              sortype = "symmetric";
    100:     else if (sym & SOR_BACKWARD_SWEEP)       sortype = "backward";
    101:     else if (sym & SOR_FORWARD_SWEEP)        sortype = "forward";
    102:     else if ((sym & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP)
    103:                                              sortype = "local_symmetric";
    104:     else if (sym & SOR_LOCAL_FORWARD_SWEEP)  sortype = "local_forward";
    105:     else if (sym & SOR_LOCAL_BACKWARD_SWEEP) sortype = "local_backward";
    106:     else                                     sortype = "unknown";
    107:     PetscViewerASCIIPrintf(viewer,"  SOR: type = %s, iterations = %D, omega = %g\n",sortype,jac->its,jac->omega);
    108:   } else {
    109:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCSOR",((PetscObject)viewer)->type_name);
    110:   }
    111:   return(0);
    112: }

    115: /* ------------------------------------------------------------------------------*/
    119: PetscErrorCode PETSCKSP_DLLEXPORT PCSORSetSymmetric_SOR(PC pc,MatSORType flag)
    120: {
    121:   PC_SOR *jac;

    124:   jac = (PC_SOR*)pc->data;
    125:   jac->sym = flag;
    126:   return(0);
    127: }

    133: PetscErrorCode PETSCKSP_DLLEXPORT PCSORSetOmega_SOR(PC pc,PetscReal omega)
    134: {
    135:   PC_SOR *jac;

    138:   if (omega >= 2.0 || omega <= 0.0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Relaxation out of range");
    139:   jac        = (PC_SOR*)pc->data;
    140:   jac->omega = omega;
    141:   return(0);
    142: }

    148: PetscErrorCode PETSCKSP_DLLEXPORT PCSORSetIterations_SOR(PC pc,PetscInt its,PetscInt lits)
    149: {
    150:   PC_SOR *jac;

    153:   jac      = (PC_SOR*)pc->data;
    154:   jac->its = its;
    155:   jac->lits = lits;
    156:   return(0);
    157: }

    160: /* ------------------------------------------------------------------------------*/
    163: /*@
    164:    PCSORSetSymmetric - Sets the SOR preconditioner to use symmetric (SSOR),
    165:    backward, or forward relaxation.  The local variants perform SOR on
    166:    each processor.  By default forward relaxation is used.

    168:    Collective on PC

    170:    Input Parameters:
    171: +  pc - the preconditioner context
    172: -  flag - one of the following
    173: .vb
    174:     SOR_FORWARD_SWEEP
    175:     SOR_BACKWARD_SWEEP
    176:     SOR_SYMMETRIC_SWEEP
    177:     SOR_LOCAL_FORWARD_SWEEP
    178:     SOR_LOCAL_BACKWARD_SWEEP
    179:     SOR_LOCAL_SYMMETRIC_SWEEP
    180: .ve

    182:    Options Database Keys:
    183: +  -pc_sor_symmetric - Activates symmetric version
    184: .  -pc_sor_backward - Activates backward version
    185: .  -pc_sor_local_forward - Activates local forward version
    186: .  -pc_sor_local_symmetric - Activates local symmetric version
    187: -  -pc_sor_local_backward - Activates local backward version

    189:    Notes:
    190:    To use the Eisenstat trick with SSOR, employ the PCEISENSTAT preconditioner,
    191:    which can be chosen with the option
    192: .  -pc_type eisenstat - Activates Eisenstat trick

    194:    Level: intermediate

    196: .keywords: PC, SOR, SSOR, set, relaxation, sweep, forward, backward, symmetric

    198: .seealso: PCEisenstatSetOmega(), PCSORSetIterations(), PCSORSetOmega()
    199: @*/
    200: PetscErrorCode PETSCKSP_DLLEXPORT PCSORSetSymmetric(PC pc,MatSORType flag)
    201: {
    202:   PetscErrorCode ierr,(*f)(PC,MatSORType);

    206:   PetscObjectQueryFunction((PetscObject)pc,"PCSORSetSymmetric_C",(void (**)(void))&f);
    207:   if (f) {
    208:     (*f)(pc,flag);
    209:   }
    210:   return(0);
    211: }

    215: /*@
    216:    PCSORSetOmega - Sets the SOR relaxation coefficient, omega
    217:    (where omega = 1.0 by default).

    219:    Collective on PC

    221:    Input Parameters:
    222: +  pc - the preconditioner context
    223: -  omega - relaxation coefficient (0 < omega < 2).

    225:    Options Database Key:
    226: .  -pc_sor_omega <omega> - Sets omega

    228:    Level: intermediate

    230: .keywords: PC, SOR, SSOR, set, relaxation, omega

    232: .seealso: PCSORSetSymmetric(), PCSORSetIterations(), PCEisenstatSetOmega()
    233: @*/
    234: PetscErrorCode PETSCKSP_DLLEXPORT PCSORSetOmega(PC pc,PetscReal omega)
    235: {
    236:   PetscErrorCode ierr,(*f)(PC,PetscReal);

    239:   PetscObjectQueryFunction((PetscObject)pc,"PCSORSetOmega_C",(void (**)(void))&f);
    240:   if (f) {
    241:     (*f)(pc,omega);
    242:   }
    243:   return(0);
    244: }

    248: /*@
    249:    PCSORSetIterations - Sets the number of inner iterations to
    250:    be used by the SOR preconditioner. The default is 1.

    252:    Collective on PC

    254:    Input Parameters:
    255: +  pc - the preconditioner context
    256: .  lits - number of local iterations, smoothings over just variables on processor
    257: -  its - number of parallel iterations to use; each parallel iteration has lits local iterations

    259:    Options Database Key:
    260: +  -pc_sor_its <its> - Sets number of iterations
    261: -  -pc_sor_lits <lits> - Sets number of local iterations

    263:    Level: intermediate

    265:    Notes: When run on one processor the number of smoothings is lits*its

    267: .keywords: PC, SOR, SSOR, set, iterations

    269: .seealso: PCSORSetOmega(), PCSORSetSymmetric()
    270: @*/
    271: PetscErrorCode PETSCKSP_DLLEXPORT PCSORSetIterations(PC pc,PetscInt its,PetscInt lits)
    272: {
    273:   PetscErrorCode ierr,(*f)(PC,PetscInt,PetscInt);

    277:   PetscObjectQueryFunction((PetscObject)pc,"PCSORSetIterations_C",(void (**)(void))&f);
    278:   if (f) {
    279:     (*f)(pc,its,lits);
    280:   }
    281:   return(0);
    282: }

    284: /*MC
    285:      PCSOR - (S)SOR (successive over relaxation, Gauss-Seidel) preconditioning

    287:    Options Database Keys:
    288: +  -pc_sor_symmetric - Activates symmetric version
    289: .  -pc_sor_backward - Activates backward version
    290: .  -pc_sor_local_forward - Activates local forward version
    291: .  -pc_sor_local_symmetric - Activates local symmetric version
    292: .  -pc_sor_local_backward - Activates local backward version
    293: .  -pc_sor_omega <omega> - Sets omega
    294: .  -pc_sor_its <its> - Sets number of iterations
    295: -  -pc_sor_lits <lits> - Sets number of local iterations

    297:    Level: beginner

    299:   Concepts: SOR, preconditioners, Gauss-Seidel

    301:    Notes: Only implemented for the AIJ  and SeqBAIJ matrix formats.
    302:           Not a true parallel SOR, in parallel this implementation corresponds to block
    303:           Jacobi with SOR on each block.

    305:           For SeqBAIJ matrices this implements point-block SOR, but the omega, its, lits options are not supported.

    307: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
    308:            PCSORSetIterations(), PCSORSetSymmetric(), PCSORSetOmega(), PCEISENSTAT
    309: M*/

    314: PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_SOR(PC pc)
    315: {
    317:   PC_SOR         *jac;

    320:   PetscNew(PC_SOR,&jac);
    321:   PetscLogObjectMemory(pc,sizeof(PC_SOR));

    323:   pc->ops->apply           = PCApply_SOR;
    324:   pc->ops->applyrichardson = PCApplyRichardson_SOR;
    325:   pc->ops->setfromoptions  = PCSetFromOptions_SOR;
    326:   pc->ops->setup           = 0;
    327:   pc->ops->view            = PCView_SOR;
    328:   pc->ops->destroy         = PCDestroy_SOR;
    329:   pc->data           = (void*)jac;
    330:   jac->sym           = SOR_LOCAL_SYMMETRIC_SWEEP;
    331:   jac->omega         = 1.0;
    332:   jac->its           = 1;
    333:   jac->lits          = 1;

    335:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCSORSetSymmetric_C","PCSORSetSymmetric_SOR",
    336:                     PCSORSetSymmetric_SOR);
    337:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCSORSetOmega_C","PCSORSetOmega_SOR",
    338:                     PCSORSetOmega_SOR);
    339:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCSORSetIterations_C","PCSORSetIterations_SOR",
    340:                     PCSORSetIterations_SOR);

    342:   return(0);
    343: }