Actual source code: vpbjacobi.c

petsc-3.10.5 2019-03-28
Report Typos and Errors

  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:   Concepts: variable point block Jacobi

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

172: M*/

174: PETSC_EXTERN PetscErrorCode PCCreate_VPBJacobi(PC pc)
175: {
176:   PC_VPBJacobi   *jac;

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

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

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