Actual source code: pbjacobi.c

      1: #define PETSCKSP_DLL

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

     8:  #include src/ksp/pc/pcimpl.h

     10: /*
     11:    Private context (data structure) for the PBJacobi preconditioner.
     12: */
     13: typedef struct {
     14:   PetscScalar *diag;
     15:   PetscInt    bs,mbs;
     16: } PC_PBJacobi;

     18: /*
     19:    Currently only implemented for baij matrices and directly access baij
     20:   data structures.
     21: */
     22:  #include src/mat/impls/baij/mpi/mpibaij.h
     23:  #include src/inline/ilu.h

     27: static PetscErrorCode PCApply_PBJacobi_2(PC pc,Vec x,Vec y)
     28: {
     29:   PC_PBJacobi    *jac = (PC_PBJacobi*)pc->data;
     31:   PetscInt       i,m = jac->mbs;
     32:   PetscScalar    *diag = jac->diag,x0,x1,*xx,*yy;
     33:
     35:   VecGetArray(x,&xx);
     36:   VecGetArray(y,&yy);
     37:   for (i=0; i<m; i++) {
     38:     x0 = xx[2*i]; x1 = xx[2*i+1];
     39:     yy[2*i]   = diag[0]*x0 + diag[2]*x1;
     40:     yy[2*i+1] = diag[1]*x0 + diag[3]*x1;
     41:     diag     += 4;
     42:   }
     43:   VecRestoreArray(x,&xx);
     44:   VecRestoreArray(y,&yy);
     45:   PetscLogFlops(6*m);
     46:   return(0);
     47: }
     50: static PetscErrorCode PCApply_PBJacobi_3(PC pc,Vec x,Vec y)
     51: {
     52:   PC_PBJacobi    *jac = (PC_PBJacobi*)pc->data;
     54:   PetscInt       i,m = jac->mbs;
     55:   PetscScalar    *diag = jac->diag,x0,x1,x2,*xx,*yy;
     56:
     58:   VecGetArray(x,&xx);
     59:   VecGetArray(y,&yy);
     60:   for (i=0; i<m; i++) {
     61:     x0 = xx[3*i]; x1 = xx[3*i+1]; x2 = xx[3*i+2];
     62:     yy[3*i]   = diag[0]*x0 + diag[3]*x1 + diag[6]*x2;
     63:     yy[3*i+1] = diag[1]*x0 + diag[4]*x1 + diag[7]*x2;
     64:     yy[3*i+2] = diag[2]*x0 + diag[5]*x1 + diag[8]*x2;
     65:     diag     += 9;
     66:   }
     67:   VecRestoreArray(x,&xx);
     68:   VecRestoreArray(y,&yy);
     69:   PetscLogFlops(15*m);
     70:   return(0);
     71: }
     74: static PetscErrorCode PCApply_PBJacobi_4(PC pc,Vec x,Vec y)
     75: {
     76:   PC_PBJacobi    *jac = (PC_PBJacobi*)pc->data;
     78:   PetscInt       i,m = jac->mbs;
     79:   PetscScalar    *diag = jac->diag,x0,x1,x2,x3,*xx,*yy;
     80:
     82:   VecGetArray(x,&xx);
     83:   VecGetArray(y,&yy);
     84:   for (i=0; i<m; i++) {
     85:     x0 = xx[4*i]; x1 = xx[4*i+1]; x2 = xx[4*i+2]; x3 = xx[4*i+3];
     86:     yy[4*i]   = diag[0]*x0 + diag[4]*x1 + diag[8]*x2  + diag[12]*x3;
     87:     yy[4*i+1] = diag[1]*x0 + diag[5]*x1 + diag[9]*x2  + diag[13]*x3;
     88:     yy[4*i+2] = diag[2]*x0 + diag[6]*x1 + diag[10]*x2 + diag[14]*x3;
     89:     yy[4*i+3] = diag[3]*x0 + diag[7]*x1 + diag[11]*x2 + diag[15]*x3;
     90:     diag     += 16;
     91:   }
     92:   VecRestoreArray(x,&xx);
     93:   VecRestoreArray(y,&yy);
     94:   PetscLogFlops(28*m);
     95:   return(0);
     96: }
     99: static PetscErrorCode PCApply_PBJacobi_5(PC pc,Vec x,Vec y)
    100: {
    101:   PC_PBJacobi    *jac = (PC_PBJacobi*)pc->data;
    103:   PetscInt       i,m = jac->mbs;
    104:   PetscScalar    *diag = jac->diag,x0,x1,x2,x3,x4,*xx,*yy;
    105:
    107:   VecGetArray(x,&xx);
    108:   VecGetArray(y,&yy);
    109:   for (i=0; i<m; i++) {
    110:     x0 = xx[5*i]; x1 = xx[5*i+1]; x2 = xx[5*i+2]; x3 = xx[5*i+3]; x4 = xx[5*i+4];
    111:     yy[5*i]   = diag[0]*x0 + diag[5]*x1 + diag[10]*x2  + diag[15]*x3 + diag[20]*x4;
    112:     yy[5*i+1] = diag[1]*x0 + diag[6]*x1 + diag[11]*x2  + diag[16]*x3 + diag[21]*x4;
    113:     yy[5*i+2] = diag[2]*x0 + diag[7]*x1 + diag[12]*x2 + diag[17]*x3 + diag[22]*x4;
    114:     yy[5*i+3] = diag[3]*x0 + diag[8]*x1 + diag[13]*x2 + diag[18]*x3 + diag[23]*x4;
    115:     yy[5*i+4] = diag[4]*x0 + diag[9]*x1 + diag[14]*x2 + diag[19]*x3 + diag[24]*x4;
    116:     diag     += 25;
    117:   }
    118:   VecRestoreArray(x,&xx);
    119:   VecRestoreArray(y,&yy);
    120:   PetscLogFlops(45*m);
    121:   return(0);
    122: }
    123: /* -------------------------------------------------------------------------- */
    126: static PetscErrorCode PCSetUp_PBJacobi(PC pc)
    127: {
    128:   PC_PBJacobi    *jac = (PC_PBJacobi*)pc->data;
    130:   PetscMPIInt    size;
    131:   PetscTruth     seqbaij,mpibaij,baij;
    132:   Mat            A = pc->pmat;
    133:   Mat_SeqBAIJ    *a;

    136:   PetscTypeCompare((PetscObject)pc->pmat,MATSEQBAIJ,&seqbaij);
    137:   PetscTypeCompare((PetscObject)pc->pmat,MATMPIBAIJ,&mpibaij);
    138:   PetscTypeCompare((PetscObject)pc->pmat,MATBAIJ,&baij);
    139:   if (!seqbaij && !mpibaij && !baij) {
    140:     SETERRQ(PETSC_ERR_SUP,"Currently only supports BAIJ matrices");
    141:   }
    142:   MPI_Comm_size(pc->comm,&size);
    143:   if (mpibaij || (baij && (size > 1))) A = ((Mat_MPIBAIJ*)A->data)->A;
    144:   if (A->m != A->n) SETERRQ(PETSC_ERR_SUP,"Supported only for square matrices and square storage");

    146:    MatSeqBAIJInvertBlockDiagonal(A);
    147:   a           = (Mat_SeqBAIJ*)A->data;
    148:   jac->diag   = a->idiag;
    149:   jac->bs     = A->bs;
    150:   jac->mbs    = a->mbs;
    151:   switch (jac->bs){
    152:     case 2:
    153:       pc->ops->apply = PCApply_PBJacobi_2;
    154:       break;
    155:     case 3:
    156:       pc->ops->apply = PCApply_PBJacobi_3;
    157:       break;
    158:     case 4:
    159:       pc->ops->apply = PCApply_PBJacobi_4;
    160:       break;
    161:     case 5:
    162:       pc->ops->apply = PCApply_PBJacobi_5;
    163:       break;
    164:     default:
    165:       SETERRQ1(PETSC_ERR_SUP,"not supported for block size %D",jac->bs);
    166:   }

    168:   return(0);
    169: }
    170: /* -------------------------------------------------------------------------- */
    173: static PetscErrorCode PCDestroy_PBJacobi(PC pc)
    174: {
    175:   PC_PBJacobi    *jac = (PC_PBJacobi*)pc->data;

    179:   /*
    180:       Free the private data structure that was hanging off the PC
    181:   */
    182:   PetscFree(jac);
    183:   return(0);
    184: }
    185: /* -------------------------------------------------------------------------- */
    186: /*MC
    187:      PCPBJACOBI - Point block Jacobi

    189:    Level: beginner

    191:   Concepts: point block Jacobi

    193:    Notes: Only implemented for the BAIJ matrix formats.

    195: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC

    197: M*/

    202: PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_PBJacobi(PC pc)
    203: {
    204:   PC_PBJacobi    *jac;

    209:   /*
    210:      Creates the private data structure for this preconditioner and
    211:      attach it to the PC object.
    212:   */
    213:   PetscNew(PC_PBJacobi,&jac);
    214:   pc->data  = (void*)jac;

    216:   /*
    217:      Logs the memory usage; this is not needed but allows PETSc to
    218:      monitor how much memory is being used for various purposes.
    219:   */
    220:   PetscLogObjectMemory(pc,sizeof(PC_PBJacobi));

    222:   /*
    223:      Initialize the pointers to vectors to ZERO; these will be used to store
    224:      diagonal entries of the matrix for fast preconditioner application.
    225:   */
    226:   jac->diag          = 0;

    228:   /*
    229:       Set the pointers for the functions that are provided above.
    230:       Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
    231:       are called, they will automatically call these functions.  Note we
    232:       choose not to provide a couple of these functions since they are
    233:       not needed.
    234:   */
    235:   pc->ops->apply               = 0; /*set depending on the block size */
    236:   pc->ops->applytranspose      = 0;
    237:   pc->ops->setup               = PCSetUp_PBJacobi;
    238:   pc->ops->destroy             = PCDestroy_PBJacobi;
    239:   pc->ops->setfromoptions      = 0;
    240:   pc->ops->view                = 0;
    241:   pc->ops->applyrichardson     = 0;
    242:   pc->ops->applysymmetricleft  = 0;
    243:   pc->ops->applysymmetricright = 0;
    244:   return(0);
    245: }