Actual source code: eisen.c

      1: #define PETSCKSP_DLL

      3: /*
      4:    Defines a  Eisenstat trick SSOR  preconditioner. This uses about
      5:  %50 of the usual amount of floating point ops used for SSOR + Krylov
      6:  method. But it requires actually solving the preconditioned problem
      7:  with both left and right preconditioning.
      8: */
     9:  #include src/ksp/pc/pcimpl.h

     11: typedef struct {
     12:   Mat        shell,A;
     13:   Vec        b,diag;     /* temporary storage for true right hand side */
     14:   PetscReal  omega;
     15:   PetscTruth usediag;    /* indicates preconditioner should include diagonal scaling*/
     16: } PC_Eisenstat;

     21: static PetscErrorCode PCMult_Eisenstat(Mat mat,Vec b,Vec x)
     22: {
     24:   PC             pc;
     25:   PC_Eisenstat   *eis;

     28:   MatShellGetContext(mat,(void **)&pc);
     29:   eis = (PC_Eisenstat*)pc->data;
     30:   MatRelax(eis->A,b,eis->omega,SOR_EISENSTAT,0.0,1,1,x);
     31:   return(0);
     32: }

     36: static PetscErrorCode PCApply_Eisenstat(PC pc,Vec x,Vec y)
     37: {
     38:   PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;

     42:   if (eis->usediag)  {VecPointwiseMult(y,x,eis->diag);}
     43:   else               {VecCopy(x,y);}
     44:   return(0);
     45: }

     49: static PetscErrorCode PCPreSolve_Eisenstat(PC pc,KSP ksp,Vec b,Vec x)
     50: {
     51:   PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;
     52:   PetscTruth     nonzero;

     56:   if (pc->mat != pc->pmat) SETERRQ(PETSC_ERR_SUP,"Cannot have different mat and pmat");
     57:
     58:   /* swap shell matrix and true matrix */
     59:   eis->A    = pc->mat;
     60:   pc->mat   = eis->shell;

     62:   if (!eis->b) {
     63:     VecDuplicate(b,&eis->b);
     64:     PetscLogObjectParent(pc,eis->b);
     65:   }
     66:
     67:   /* save true b, other option is to swap pointers */
     68:   VecCopy(b,eis->b);

     70:   /* if nonzero initial guess, modify x */
     71:   KSPGetInitialGuessNonzero(ksp,&nonzero);
     72:   if (nonzero) {
     73:     MatRelax(eis->A,x,eis->omega,SOR_APPLY_UPPER,0.0,1,1,x);
     74:   }

     76:   /* modify b by (L + D)^{-1} */
     77:     MatRelax(eis->A,b,eis->omega,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_FORWARD_SWEEP),0.0,1,1,b);
     78:   return(0);
     79: }

     83: static PetscErrorCode PCPostSolve_Eisenstat(PC pc,KSP ksp,Vec b,Vec x)
     84: {
     85:   PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;

     89:     MatRelax(eis->A,x,eis->omega,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_BACKWARD_SWEEP),0.0,1,1,x);
     90:   pc->mat = eis->A;
     91:   /* get back true b */
     92:   VecCopy(eis->b,b);
     93:   return(0);
     94: }

     98: static PetscErrorCode PCDestroy_Eisenstat(PC pc)
     99: {
    100:   PC_Eisenstat   *eis = (PC_Eisenstat *)pc->data;

    104:   if (eis->b)     {VecDestroy(eis->b);}
    105:   if (eis->shell) {MatDestroy(eis->shell);}
    106:   if (eis->diag)  {VecDestroy(eis->diag);}
    107:   PetscFree(eis);
    108:   return(0);
    109: }

    113: static PetscErrorCode PCSetFromOptions_Eisenstat(PC pc)
    114: {
    115:   PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;
    117:   PetscTruth     flg;

    120:   PetscOptionsHead("Eisenstat SSOR options");
    121:     PetscOptionsReal("-pc_eisenstat_omega","Relaxation factor 0 < omega < 2","PCEisenstatSetOmega",eis->omega,&eis->omega,0);
    122:     PetscOptionsName("-pc_eisenstat_no_diagonal_scaling","Do not use standard diagonal scaling","PCEisenstatNoDiagonalScaling",&flg);
    123:     if (flg) {
    124:       PCEisenstatNoDiagonalScaling(pc);
    125:     }
    126:   PetscOptionsTail();
    127:   return(0);
    128: }

    132: static PetscErrorCode PCView_Eisenstat(PC pc,PetscViewer viewer)
    133: {
    134:   PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;
    136:   PetscTruth     iascii;

    139:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
    140:   if (iascii) {
    141:     PetscViewerASCIIPrintf(viewer,"Eisenstat: omega = %g\n",eis->omega);
    142:     if (eis->usediag) {
    143:       PetscViewerASCIIPrintf(viewer,"Eisenstat: Using diagonal scaling (default)\n");
    144:     } else {
    145:       PetscViewerASCIIPrintf(viewer,"Eisenstat: Not using diagonal scaling\n");
    146:     }
    147:   } else {
    148:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for Eisenstat PC",((PetscObject)viewer)->type_name);
    149:   }
    150:   return(0);
    151: }

    155: static PetscErrorCode PCSetUp_Eisenstat(PC pc)
    156: {
    158:   PetscInt       M,N,m,n;
    159:   PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;

    162:   if (!pc->setupcalled) {
    163:     MatGetSize(pc->mat,&M,&N);
    164:     MatGetLocalSize(pc->mat,&m,&n);
    165:     MatCreate(pc->comm,&eis->shell);
    166:     MatSetSizes(eis->shell,m,N,M,N);
    167:     MatSetType(eis->shell,MATSHELL);
    168:     MatShellSetContext(eis->shell,(void*)pc);
    169:     PetscLogObjectParent(pc,eis->shell);
    170:     MatShellSetOperation(eis->shell,MATOP_MULT,(void(*)(void))PCMult_Eisenstat);
    171:   }
    172:   if (!eis->usediag) return(0);
    173:   if (!pc->setupcalled) {
    174:     MatGetVecs(pc->pmat,&eis->diag,0);
    175:     PetscLogObjectParent(pc,eis->diag);
    176:   }
    177:   MatGetDiagonal(pc->pmat,eis->diag);
    178:   return(0);
    179: }

    181: /* --------------------------------------------------------------------*/

    186: PetscErrorCode PETSCKSP_DLLEXPORT PCEisenstatSetOmega_Eisenstat(PC pc,PetscReal omega)
    187: {
    188:   PC_Eisenstat  *eis;

    191:   if (omega >= 2.0 || omega <= 0.0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Relaxation out of range");
    192:   eis = (PC_Eisenstat*)pc->data;
    193:   eis->omega = omega;
    194:   return(0);
    195: }

    201: PetscErrorCode PETSCKSP_DLLEXPORT PCEisenstatNoDiagonalScaling_Eisenstat(PC pc)
    202: {
    203:   PC_Eisenstat *eis;

    206:   eis = (PC_Eisenstat*)pc->data;
    207:   eis->usediag = PETSC_FALSE;
    208:   return(0);
    209: }

    214: /*@
    215:    PCEisenstatSetOmega - Sets the SSOR relaxation coefficient, omega,
    216:    to use with Eisenstat's trick (where omega = 1.0 by default).

    218:    Collective on PC

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

    224:    Options Database Key:
    225: .  -pc_eisenstat_omega <omega> - Sets omega

    227:    Notes:
    228:    The Eisenstat trick implementation of SSOR requires about 50% of the
    229:    usual amount of floating point operations used for SSOR + Krylov method;
    230:    however, the preconditioned problem must be solved with both left
    231:    and right preconditioning.

    233:    To use SSOR without the Eisenstat trick, employ the PCSOR preconditioner,
    234:    which can be chosen with the database options
    235: $    -pc_type  sor  -pc_sor_symmetric

    237:    Level: intermediate

    239: .keywords: PC, Eisenstat, set, SOR, SSOR, relaxation, omega

    241: .seealso: PCSORSetOmega()
    242: @*/
    243: PetscErrorCode PETSCKSP_DLLEXPORT PCEisenstatSetOmega(PC pc,PetscReal omega)
    244: {
    245:   PetscErrorCode ierr,(*f)(PC,PetscReal);

    249:   PetscObjectQueryFunction((PetscObject)pc,"PCEisenstatSetOmega_C",(void (**)(void))&f);
    250:   if (f) {
    251:     (*f)(pc,omega);
    252:   }
    253:   return(0);
    254: }

    258: /*@
    259:    PCEisenstatNoDiagonalScaling - Causes the Eisenstat preconditioner
    260:    not to do additional diagonal preconditioning. For matrices with a constant
    261:    along the diagonal, this may save a small amount of work.

    263:    Collective on PC

    265:    Input Parameter:
    266: .  pc - the preconditioner context

    268:    Options Database Key:
    269: .  -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatNoDiagonalScaling()

    271:    Level: intermediate

    273:    Note:
    274:      If you use the KPSSetDiagonalScaling() or -ksp_diagonal_scale option then you will
    275:    likley want to use this routine since it will save you some unneeded flops.

    277: .keywords: PC, Eisenstat, use, diagonal, scaling, SSOR

    279: .seealso: PCEisenstatSetOmega()
    280: @*/
    281: PetscErrorCode PETSCKSP_DLLEXPORT PCEisenstatNoDiagonalScaling(PC pc)
    282: {
    283:   PetscErrorCode ierr,(*f)(PC);

    287:   PetscObjectQueryFunction((PetscObject)pc,"PCEisenstatNoDiagonalScaling_C",(void (**)(void))&f);
    288:   if (f) {
    289:     (*f)(pc);
    290:   }
    291:   return(0);
    292: }

    294: /* --------------------------------------------------------------------*/

    296: /*MC
    297:      PCEISENSTAT - An implementation of SSOR (symmetric successive over relaxation, symmetric Gauss-Seidel)
    298:            preconditioning that incorporates Eisenstat's trick to reduce the amount of computation needed.

    300:    Options Database Keys:
    301: +  -pc_eisenstat_omega <omega> - Sets omega
    302: -  -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatNoDiagonalScaling()

    304:    Level: beginner

    306:   Concepts: SOR, preconditioners, Gauss-Seidel, Eisenstat's trick

    308:    Notes: Only implemented for the SeqAIJ matrix format.
    309:           Not a true parallel SOR, in parallel this implementation corresponds to block
    310:           Jacobi with SOR on each block.

    312: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
    313:            PCEisenstatNoDiagonalScaling(), PCEisenstatSetOmega(), PCSOR
    314: M*/

    319: PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_Eisenstat(PC pc)
    320: {
    322:   PC_Eisenstat   *eis;

    325:   PetscNew(PC_Eisenstat,&eis);
    326:   PetscLogObjectMemory(pc,sizeof(PC_Eisenstat));

    328:   pc->ops->apply           = PCApply_Eisenstat;
    329:   pc->ops->presolve        = PCPreSolve_Eisenstat;
    330:   pc->ops->postsolve       = PCPostSolve_Eisenstat;
    331:   pc->ops->applyrichardson = 0;
    332:   pc->ops->setfromoptions  = PCSetFromOptions_Eisenstat;
    333:   pc->ops->destroy         = PCDestroy_Eisenstat;
    334:   pc->ops->view            = PCView_Eisenstat;
    335:   pc->ops->setup           = PCSetUp_Eisenstat;

    337:   pc->data           = (void*)eis;
    338:   eis->omega         = 1.0;
    339:   eis->b             = 0;
    340:   eis->diag          = 0;
    341:   eis->usediag       = PETSC_TRUE;

    343:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCEisenstatSetOmega_C","PCEisenstatSetOmega_Eisenstat",
    344:                     PCEisenstatSetOmega_Eisenstat);
    345:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCEisenstatNoDiagonalScaling_C",
    346:                     "PCEisenstatNoDiagonalScaling_Eisenstat",
    347:                     PCEisenstatNoDiagonalScaling_Eisenstat);
    348:  return(0);
    349: }