Actual source code: vpbjacobi.c

      2: /*
      3:    Include files needed for the variable size block PBJacobi preconditioner:
      4:      pcimpl.h - private include file intended for use by all preconditioners
      5: */

      7: #include <petsc/private/pcimpl.h>

      9: /*
     10:    Private context (data structure) for the VPBJacobi preconditioner.
     11: */
     12: typedef struct {
     13:   MatScalar *diag;
     14: } PC_VPBJacobi;

     17: static PetscErrorCode PCApply_VPBJacobi(PC pc,Vec x,Vec y)
     18: {
     19:   PC_VPBJacobi      *jac = (PC_VPBJacobi*)pc->data;
     20:   PetscErrorCode    ierr;
     21:   PetscInt          i,ncnt = 0;
     22:   const MatScalar   *diag = jac->diag;
     23:   PetscInt          ib,jb,bs;
     24:   const PetscScalar *xx;
     25:   PetscScalar       *yy,x0,x1,x2,x3,x4,x5,x6;
     26:   PetscInt          nblocks;
     27:   const PetscInt    *bsizes;

     30:   MatGetVariableBlockSizes(pc->pmat,&nblocks,&bsizes);
     31:   VecGetArrayRead(x,&xx);
     32:   VecGetArray(y,&yy);
     33:   for (i=0; i<nblocks; i++) {
     34:     bs = bsizes[i];
     35:     switch (bs) {
     36:     case 1:
     37:       yy[ncnt] = *diag*xx[ncnt];
     38:       break;
     39:     case 2:
     40:       x0         = xx[ncnt]; x1 = xx[ncnt+1];
     41:       yy[ncnt]   = diag[0]*x0 + diag[2]*x1;
     42:       yy[ncnt+1] = diag[1]*x0 + diag[3]*x1;
     43:       break;
     44:     case 3:
     45:       x0 = xx[ncnt]; x1 = xx[ncnt+1]; x2 = xx[ncnt+2];
     46:       yy[ncnt]   = diag[0]*x0 + diag[3]*x1 + diag[6]*x2;
     47:       yy[ncnt+1] = diag[1]*x0 + diag[4]*x1 + diag[7]*x2;
     48:       yy[ncnt+2] = diag[2]*x0 + diag[5]*x1 + diag[8]*x2;
     49:       break;
     50:     case 4:
     51:       x0 = xx[ncnt]; x1 = xx[ncnt+1]; x2 = xx[ncnt+2]; x3 = xx[ncnt+3];
     52:       yy[ncnt]   = diag[0]*x0 + diag[4]*x1 + diag[8]*x2  + diag[12]*x3;
     53:       yy[ncnt+1] = diag[1]*x0 + diag[5]*x1 + diag[9]*x2  + diag[13]*x3;
     54:       yy[ncnt+2] = diag[2]*x0 + diag[6]*x1 + diag[10]*x2 + diag[14]*x3;
     55:       yy[ncnt+3] = diag[3]*x0 + diag[7]*x1 + diag[11]*x2 + diag[15]*x3;
     56:       break;
     57:     case 5:
     58:       x0 = xx[ncnt]; x1 = xx[ncnt+1]; x2 = xx[ncnt+2]; x3 = xx[ncnt+3]; x4 = xx[ncnt+4];
     59:       yy[ncnt]   = diag[0]*x0 + diag[5]*x1 + diag[10]*x2  + diag[15]*x3 + diag[20]*x4;
     60:       yy[ncnt+1] = diag[1]*x0 + diag[6]*x1 + diag[11]*x2  + diag[16]*x3 + diag[21]*x4;
     61:       yy[ncnt+2] = diag[2]*x0 + diag[7]*x1 + diag[12]*x2 + diag[17]*x3 + diag[22]*x4;
     62:       yy[ncnt+3] = diag[3]*x0 + diag[8]*x1 + diag[13]*x2 + diag[18]*x3 + diag[23]*x4;
     63:       yy[ncnt+4] = diag[4]*x0 + diag[9]*x1 + diag[14]*x2 + diag[19]*x3 + diag[24]*x4;
     64:       break;
     65:     case 6:
     66:       x0 = xx[ncnt]; x1 = xx[ncnt+1]; x2 = xx[ncnt+2]; x3 = xx[ncnt+3]; x4 = xx[ncnt+4]; x5 = xx[ncnt+5];
     67:       yy[ncnt]   = diag[0]*x0 + diag[6]*x1  + diag[12]*x2  + diag[18]*x3 + diag[24]*x4 + diag[30]*x5;
     68:       yy[ncnt+1] = diag[1]*x0 + diag[7]*x1  + diag[13]*x2  + diag[19]*x3 + diag[25]*x4 + diag[31]*x5;
     69:       yy[ncnt+2] = diag[2]*x0 + diag[8]*x1  + diag[14]*x2  + diag[20]*x3 + diag[26]*x4 + diag[32]*x5;
     70:       yy[ncnt+3] = diag[3]*x0 + diag[9]*x1  + diag[15]*x2  + diag[21]*x3 + diag[27]*x4 + diag[33]*x5;
     71:       yy[ncnt+4] = diag[4]*x0 + diag[10]*x1 + diag[16]*x2  + diag[22]*x3 + diag[28]*x4 + diag[34]*x5;
     72:       yy[ncnt+5] = diag[5]*x0 + diag[11]*x1 + diag[17]*x2  + diag[23]*x3 + diag[29]*x4 + diag[35]*x5;
     73:       break;
     74:     case 7:
     75:       x0 = xx[ncnt]; x1 = xx[ncnt+1]; x2 = xx[ncnt+2]; x3 = xx[ncnt+3]; x4 = xx[ncnt+4]; x5 = xx[ncnt+5]; x6 = xx[ncnt+6];
     76:       yy[ncnt]   = diag[0]*x0 + diag[7]*x1  + diag[14]*x2  + diag[21]*x3 + diag[28]*x4 + diag[35]*x5 + diag[42]*x6;
     77:       yy[ncnt+1] = diag[1]*x0 + diag[8]*x1  + diag[15]*x2  + diag[22]*x3 + diag[29]*x4 + diag[36]*x5 + diag[43]*x6;
     78:       yy[ncnt+2] = diag[2]*x0 + diag[9]*x1  + diag[16]*x2  + diag[23]*x3 + diag[30]*x4 + diag[37]*x5 + diag[44]*x6;
     79:       yy[ncnt+3] = diag[3]*x0 + diag[10]*x1 + diag[17]*x2  + diag[24]*x3 + diag[31]*x4 + diag[38]*x5 + diag[45]*x6;
     80:       yy[ncnt+4] = diag[4]*x0 + diag[11]*x1 + diag[18]*x2  + diag[25]*x3 + diag[32]*x4 + diag[39]*x5 + diag[46]*x6;
     81:       yy[ncnt+5] = diag[5]*x0 + diag[12]*x1 + diag[19]*x2  + diag[26]*x3 + diag[33]*x4 + diag[40]*x5 + diag[47]*x6;
     82:       yy[ncnt+6] = diag[6]*x0 + diag[13]*x1 + diag[20]*x2  + diag[27]*x3 + diag[34]*x4 + diag[41]*x5 + diag[48]*x6;
     83:       break;
     84:     default:
     85:       for (ib=0; ib<bs; ib++){
     86:         PetscScalar rowsum = 0;
     87:         for (jb=0; jb<bs; jb++){
     88:           rowsum += diag[ib+jb*bs] * xx[ncnt+jb];
     89:         }
     90:         yy[ncnt+ib] = rowsum;
     91:       }
     92:     }
     93:     ncnt += bsizes[i];
     94:     diag += bsizes[i]*bsizes[i];
     95:   }
     96:   VecRestoreArrayRead(x,&xx);
     97:   VecRestoreArray(y,&yy);
     98:   return(0);
     99: }

    103: /* -------------------------------------------------------------------------- */
    104: static PetscErrorCode PCSetUp_VPBJacobi(PC pc)
    105: {
    106:   PC_VPBJacobi    *jac = (PC_VPBJacobi*)pc->data;
    108:   Mat            A = pc->pmat;
    109:   MatFactorError err;
    110:   PetscInt       i,nsize = 0,nlocal;
    111:   PetscInt       nblocks;
    112:   const PetscInt *bsizes;

    115:   MatGetVariableBlockSizes(pc->pmat,&nblocks,&bsizes);
    116:   MatGetLocalSize(pc->pmat,&nlocal,NULL);
    117:   if (nlocal && !nblocks) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call MatSetVariableBlockSizes() before using PCVPBJACOBI");
    118:   if (!jac->diag) {
    119:     for (i=0; i<nblocks; i++) nsize += bsizes[i]*bsizes[i];
    120:     PetscMalloc1(nsize,&jac->diag);
    121:   }
    122:   MatInvertVariableBlockDiagonal(A,nblocks,bsizes,jac->diag);
    123:   MatFactorGetError(A,&err);
    124:   if (err) pc->failedreason = (PCFailedReason)err;
    125:   pc->ops->apply = PCApply_VPBJacobi;
    126:   return(0);
    127: }
    128: /* -------------------------------------------------------------------------- */
    129: static PetscErrorCode PCDestroy_VPBJacobi(PC pc)
    130: {
    131:   PC_VPBJacobi    *jac = (PC_VPBJacobi*)pc->data;

    135:   /*
    136:       Free the private data structure that was hanging off the PC
    137:   */
    138:   PetscFree(jac->diag);
    139:   PetscFree(pc->data);
    140:   return(0);
    141: }

    143: /* -------------------------------------------------------------------------- */
    144: /*MC
    145:      PCVPBJACOBI - Variable size point block Jacobi preconditioner

    148:    Notes:
    149:     See PCJACOBI for point Jacobi preconditioning, PCPBJACOBI for fixed point block size, and PCBJACOBI for large size blocks

    151:    This works for AIJ matrices

    153:    Uses dense LU factorization with partial pivoting to invert the blocks; if a zero pivot
    154:    is detected a PETSc error is generated.

    156:    One must call MatSetVariableBlockSizes() to use this preconditioner
    157:    Developer Notes:
    158:     This should support the PCSetErrorIfFailure() flag set to PETSC_TRUE to allow
    159:    the factorization to continue even after a zero pivot is found resulting in a Nan and hence
    160:    terminating KSP with a KSP_DIVERGED_NANORIF allowing
    161:    a nonlinear solver/ODE integrator to recover without stopping the program as currently happens.

    163:    Perhaps should provide an option that allows generation of a valid preconditioner
    164:    even if a block is singular as the PCJACOBI does.

    166:    Level: beginner

    168: .seealso:  MatSetVariableBlockSizes(), PCCreate(), PCSetType(), PCType (for list of available types), PC, PCJACOBI

    170: M*/

    172: PETSC_EXTERN PetscErrorCode PCCreate_VPBJacobi(PC pc)
    173: {
    174:   PC_VPBJacobi   *jac;

    178:   /*
    179:      Creates the private data structure for this preconditioner and
    180:      attach it to the PC object.
    181:   */
    182:   PetscNewLog(pc,&jac);
    183:   pc->data = (void*)jac;

    185:   /*
    186:      Initialize the pointers to vectors to ZERO; these will be used to store
    187:      diagonal entries of the matrix for fast preconditioner application.
    188:   */
    189:   jac->diag = NULL;

    191:   /*
    192:       Set the pointers for the functions that are provided above.
    193:       Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
    194:       are called, they will automatically call these functions.  Note we
    195:       choose not to provide a couple of these functions since they are
    196:       not needed.
    197:   */
    198:   pc->ops->apply               = PCApply_VPBJacobi;
    199:   pc->ops->applytranspose      = NULL;
    200:   pc->ops->setup               = PCSetUp_VPBJacobi;
    201:   pc->ops->destroy             = PCDestroy_VPBJacobi;
    202:   pc->ops->setfromoptions      = NULL;
    203:   pc->ops->applyrichardson     = NULL;
    204:   pc->ops->applysymmetricleft  = NULL;
    205:   pc->ops->applysymmetricright = NULL;
    206:   return(0);
    207: }