1 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 */ 6 7 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 8 9 /* 10 Private context (data structure) for the VPBJacobi preconditioner. 11 */ 12 typedef struct { 13 MatScalar *diag; 14 } PC_VPBJacobi; 15 16 static PetscErrorCode PCApply_VPBJacobi(PC pc,Vec x,Vec y) 17 { 18 PC_VPBJacobi *jac = (PC_VPBJacobi*)pc->data; 19 PetscInt i,ncnt = 0; 20 const MatScalar *diag = jac->diag; 21 PetscInt ib,jb,bs; 22 const PetscScalar *xx; 23 PetscScalar *yy,x0,x1,x2,x3,x4,x5,x6; 24 PetscInt nblocks; 25 const PetscInt *bsizes; 26 27 PetscFunctionBegin; 28 PetscCall(MatGetVariableBlockSizes(pc->pmat,&nblocks,&bsizes)); 29 PetscCall(VecGetArrayRead(x,&xx)); 30 PetscCall(VecGetArray(y,&yy)); 31 for (i=0; i<nblocks; i++) { 32 bs = bsizes[i]; 33 switch (bs) { 34 case 1: 35 yy[ncnt] = *diag*xx[ncnt]; 36 break; 37 case 2: 38 x0 = xx[ncnt]; x1 = xx[ncnt+1]; 39 yy[ncnt] = diag[0]*x0 + diag[2]*x1; 40 yy[ncnt+1] = diag[1]*x0 + diag[3]*x1; 41 break; 42 case 3: 43 x0 = xx[ncnt]; x1 = xx[ncnt+1]; x2 = xx[ncnt+2]; 44 yy[ncnt] = diag[0]*x0 + diag[3]*x1 + diag[6]*x2; 45 yy[ncnt+1] = diag[1]*x0 + diag[4]*x1 + diag[7]*x2; 46 yy[ncnt+2] = diag[2]*x0 + diag[5]*x1 + diag[8]*x2; 47 break; 48 case 4: 49 x0 = xx[ncnt]; x1 = xx[ncnt+1]; x2 = xx[ncnt+2]; x3 = xx[ncnt+3]; 50 yy[ncnt] = diag[0]*x0 + diag[4]*x1 + diag[8]*x2 + diag[12]*x3; 51 yy[ncnt+1] = diag[1]*x0 + diag[5]*x1 + diag[9]*x2 + diag[13]*x3; 52 yy[ncnt+2] = diag[2]*x0 + diag[6]*x1 + diag[10]*x2 + diag[14]*x3; 53 yy[ncnt+3] = diag[3]*x0 + diag[7]*x1 + diag[11]*x2 + diag[15]*x3; 54 break; 55 case 5: 56 x0 = xx[ncnt]; x1 = xx[ncnt+1]; x2 = xx[ncnt+2]; x3 = xx[ncnt+3]; x4 = xx[ncnt+4]; 57 yy[ncnt] = diag[0]*x0 + diag[5]*x1 + diag[10]*x2 + diag[15]*x3 + diag[20]*x4; 58 yy[ncnt+1] = diag[1]*x0 + diag[6]*x1 + diag[11]*x2 + diag[16]*x3 + diag[21]*x4; 59 yy[ncnt+2] = diag[2]*x0 + diag[7]*x1 + diag[12]*x2 + diag[17]*x3 + diag[22]*x4; 60 yy[ncnt+3] = diag[3]*x0 + diag[8]*x1 + diag[13]*x2 + diag[18]*x3 + diag[23]*x4; 61 yy[ncnt+4] = diag[4]*x0 + diag[9]*x1 + diag[14]*x2 + diag[19]*x3 + diag[24]*x4; 62 break; 63 case 6: 64 x0 = xx[ncnt]; x1 = xx[ncnt+1]; x2 = xx[ncnt+2]; x3 = xx[ncnt+3]; x4 = xx[ncnt+4]; x5 = xx[ncnt+5]; 65 yy[ncnt] = diag[0]*x0 + diag[6]*x1 + diag[12]*x2 + diag[18]*x3 + diag[24]*x4 + diag[30]*x5; 66 yy[ncnt+1] = diag[1]*x0 + diag[7]*x1 + diag[13]*x2 + diag[19]*x3 + diag[25]*x4 + diag[31]*x5; 67 yy[ncnt+2] = diag[2]*x0 + diag[8]*x1 + diag[14]*x2 + diag[20]*x3 + diag[26]*x4 + diag[32]*x5; 68 yy[ncnt+3] = diag[3]*x0 + diag[9]*x1 + diag[15]*x2 + diag[21]*x3 + diag[27]*x4 + diag[33]*x5; 69 yy[ncnt+4] = diag[4]*x0 + diag[10]*x1 + diag[16]*x2 + diag[22]*x3 + diag[28]*x4 + diag[34]*x5; 70 yy[ncnt+5] = diag[5]*x0 + diag[11]*x1 + diag[17]*x2 + diag[23]*x3 + diag[29]*x4 + diag[35]*x5; 71 break; 72 case 7: 73 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]; 74 yy[ncnt] = diag[0]*x0 + diag[7]*x1 + diag[14]*x2 + diag[21]*x3 + diag[28]*x4 + diag[35]*x5 + diag[42]*x6; 75 yy[ncnt+1] = diag[1]*x0 + diag[8]*x1 + diag[15]*x2 + diag[22]*x3 + diag[29]*x4 + diag[36]*x5 + diag[43]*x6; 76 yy[ncnt+2] = diag[2]*x0 + diag[9]*x1 + diag[16]*x2 + diag[23]*x3 + diag[30]*x4 + diag[37]*x5 + diag[44]*x6; 77 yy[ncnt+3] = diag[3]*x0 + diag[10]*x1 + diag[17]*x2 + diag[24]*x3 + diag[31]*x4 + diag[38]*x5 + diag[45]*x6; 78 yy[ncnt+4] = diag[4]*x0 + diag[11]*x1 + diag[18]*x2 + diag[25]*x3 + diag[32]*x4 + diag[39]*x5 + diag[46]*x6; 79 yy[ncnt+5] = diag[5]*x0 + diag[12]*x1 + diag[19]*x2 + diag[26]*x3 + diag[33]*x4 + diag[40]*x5 + diag[47]*x6; 80 yy[ncnt+6] = diag[6]*x0 + diag[13]*x1 + diag[20]*x2 + diag[27]*x3 + diag[34]*x4 + diag[41]*x5 + diag[48]*x6; 81 break; 82 default: 83 for (ib=0; ib<bs; ib++) { 84 PetscScalar rowsum = 0; 85 for (jb=0; jb<bs; jb++) { 86 rowsum += diag[ib+jb*bs] * xx[ncnt+jb]; 87 } 88 yy[ncnt+ib] = rowsum; 89 } 90 } 91 ncnt += bsizes[i]; 92 diag += bsizes[i]*bsizes[i]; 93 } 94 PetscCall(VecRestoreArrayRead(x,&xx)); 95 PetscCall(VecRestoreArray(y,&yy)); 96 PetscFunctionReturn(0); 97 } 98 99 /* -------------------------------------------------------------------------- */ 100 static PetscErrorCode PCSetUp_VPBJacobi(PC pc) 101 { 102 PC_VPBJacobi *jac = (PC_VPBJacobi*)pc->data; 103 Mat A = pc->pmat; 104 MatFactorError err; 105 PetscInt i,nsize = 0,nlocal; 106 PetscInt nblocks; 107 const PetscInt *bsizes; 108 109 PetscFunctionBegin; 110 PetscCall(MatGetVariableBlockSizes(pc->pmat,&nblocks,&bsizes)); 111 PetscCall(MatGetLocalSize(pc->pmat,&nlocal,NULL)); 112 PetscCheck(!nlocal || nblocks,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call MatSetVariableBlockSizes() before using PCVPBJACOBI"); 113 if (!jac->diag) { 114 for (i=0; i<nblocks; i++) nsize += bsizes[i]*bsizes[i]; 115 PetscCall(PetscMalloc1(nsize,&jac->diag)); 116 } 117 PetscCall(MatInvertVariableBlockDiagonal(A,nblocks,bsizes,jac->diag)); 118 PetscCall(MatFactorGetError(A,&err)); 119 if (err) pc->failedreason = (PCFailedReason)err; 120 pc->ops->apply = PCApply_VPBJacobi; 121 PetscFunctionReturn(0); 122 } 123 /* -------------------------------------------------------------------------- */ 124 static PetscErrorCode PCDestroy_VPBJacobi(PC pc) 125 { 126 PC_VPBJacobi *jac = (PC_VPBJacobi*)pc->data; 127 128 PetscFunctionBegin; 129 /* 130 Free the private data structure that was hanging off the PC 131 */ 132 PetscCall(PetscFree(jac->diag)); 133 PetscCall(PetscFree(pc->data)); 134 PetscFunctionReturn(0); 135 } 136 137 /* -------------------------------------------------------------------------- */ 138 /*MC 139 PCVPBJACOBI - Variable size point block Jacobi preconditioner 140 141 Notes: 142 See PCJACOBI for point Jacobi preconditioning, PCPBJACOBI for fixed point block size, and PCBJACOBI for large size blocks 143 144 This works for AIJ matrices 145 146 Uses dense LU factorization with partial pivoting to invert the blocks; if a zero pivot 147 is detected a PETSc error is generated. 148 149 One must call MatSetVariableBlockSizes() to use this preconditioner 150 151 Developer Notes: 152 This should support the PCSetErrorIfFailure() flag set to PETSC_TRUE to allow 153 the factorization to continue even after a zero pivot is found resulting in a Nan and hence 154 terminating KSP with a KSP_DIVERGED_NANORIF allowing 155 a nonlinear solver/ODE integrator to recover without stopping the program as currently happens. 156 157 Perhaps should provide an option that allows generation of a valid preconditioner 158 even if a block is singular as the PCJACOBI does. 159 160 Level: beginner 161 162 .seealso: MatSetVariableBlockSizes(), PCCreate(), PCSetType(), PCType (for list of available types), PC, PCJACOBI, PCPBJACOBI, PCBJACOBI 163 164 M*/ 165 166 PETSC_EXTERN PetscErrorCode PCCreate_VPBJacobi(PC pc) 167 { 168 PC_VPBJacobi *jac; 169 170 PetscFunctionBegin; 171 /* 172 Creates the private data structure for this preconditioner and 173 attach it to the PC object. 174 */ 175 PetscCall(PetscNewLog(pc,&jac)); 176 pc->data = (void*)jac; 177 178 /* 179 Initialize the pointers to vectors to ZERO; these will be used to store 180 diagonal entries of the matrix for fast preconditioner application. 181 */ 182 jac->diag = NULL; 183 184 /* 185 Set the pointers for the functions that are provided above. 186 Now when the user-level routines (such as PCApply(), PCDestroy(), etc.) 187 are called, they will automatically call these functions. Note we 188 choose not to provide a couple of these functions since they are 189 not needed. 190 */ 191 pc->ops->apply = PCApply_VPBJacobi; 192 pc->ops->applytranspose = NULL; 193 pc->ops->setup = PCSetUp_VPBJacobi; 194 pc->ops->destroy = PCDestroy_VPBJacobi; 195 pc->ops->setfromoptions = NULL; 196 pc->ops->applyrichardson = NULL; 197 pc->ops->applysymmetricleft = NULL; 198 pc->ops->applysymmetricright = NULL; 199 PetscFunctionReturn(0); 200 } 201