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 static PetscErrorCode PCSetUp_VPBJacobi(PC pc) 100 { 101 PC_VPBJacobi *jac = (PC_VPBJacobi*)pc->data; 102 Mat A = pc->pmat; 103 MatFactorError err; 104 PetscInt i,nsize = 0,nlocal; 105 PetscInt nblocks; 106 const PetscInt *bsizes; 107 108 PetscFunctionBegin; 109 PetscCall(MatGetVariableBlockSizes(pc->pmat,&nblocks,&bsizes)); 110 PetscCall(MatGetLocalSize(pc->pmat,&nlocal,NULL)); 111 PetscCheck(!nlocal || nblocks,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call MatSetVariableBlockSizes() before using PCVPBJACOBI"); 112 if (!jac->diag) { 113 for (i=0; i<nblocks; i++) nsize += bsizes[i]*bsizes[i]; 114 PetscCall(PetscMalloc1(nsize,&jac->diag)); 115 } 116 PetscCall(MatInvertVariableBlockDiagonal(A,nblocks,bsizes,jac->diag)); 117 PetscCall(MatFactorGetError(A,&err)); 118 if (err) pc->failedreason = (PCFailedReason)err; 119 pc->ops->apply = PCApply_VPBJacobi; 120 PetscFunctionReturn(0); 121 } 122 123 static PetscErrorCode PCDestroy_VPBJacobi(PC pc) 124 { 125 PC_VPBJacobi *jac = (PC_VPBJacobi*)pc->data; 126 127 PetscFunctionBegin; 128 /* 129 Free the private data structure that was hanging off the PC 130 */ 131 PetscCall(PetscFree(jac->diag)); 132 PetscCall(PetscFree(pc->data)); 133 PetscFunctionReturn(0); 134 } 135 136 /*MC 137 PCVPBJACOBI - Variable size point block Jacobi preconditioner 138 139 Notes: 140 See PCJACOBI for point Jacobi preconditioning, PCPBJACOBI for fixed point block size, and PCBJACOBI for large size blocks 141 142 This works for AIJ matrices 143 144 Uses dense LU factorization with partial pivoting to invert the blocks; if a zero pivot 145 is detected a PETSc error is generated. 146 147 One must call MatSetVariableBlockSizes() to use this preconditioner 148 149 Developer Notes: 150 This should support the PCSetErrorIfFailure() flag set to PETSC_TRUE to allow 151 the factorization to continue even after a zero pivot is found resulting in a Nan and hence 152 terminating KSP with a KSP_DIVERGED_NANORIF allowing 153 a nonlinear solver/ODE integrator to recover without stopping the program as currently happens. 154 155 Perhaps should provide an option that allows generation of a valid preconditioner 156 even if a block is singular as the PCJACOBI does. 157 158 Level: beginner 159 160 .seealso: `MatSetVariableBlockSizes()`, `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCJACOBI`, `PCPBJACOBI`, `PCBJACOBI` 161 162 M*/ 163 164 PETSC_EXTERN PetscErrorCode PCCreate_VPBJacobi(PC pc) 165 { 166 PC_VPBJacobi *jac; 167 168 PetscFunctionBegin; 169 /* 170 Creates the private data structure for this preconditioner and 171 attach it to the PC object. 172 */ 173 PetscCall(PetscNewLog(pc,&jac)); 174 pc->data = (void*)jac; 175 176 /* 177 Initialize the pointers to vectors to ZERO; these will be used to store 178 diagonal entries of the matrix for fast preconditioner application. 179 */ 180 jac->diag = NULL; 181 182 /* 183 Set the pointers for the functions that are provided above. 184 Now when the user-level routines (such as PCApply(), PCDestroy(), etc.) 185 are called, they will automatically call these functions. Note we 186 choose not to provide a couple of these functions since they are 187 not needed. 188 */ 189 pc->ops->apply = PCApply_VPBJacobi; 190 pc->ops->applytranspose = NULL; 191 pc->ops->setup = PCSetUp_VPBJacobi; 192 pc->ops->destroy = PCDestroy_VPBJacobi; 193 pc->ops->setfromoptions = NULL; 194 pc->ops->applyrichardson = NULL; 195 pc->ops->applysymmetricleft = NULL; 196 pc->ops->applysymmetricright = NULL; 197 PetscFunctionReturn(0); 198 } 199