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 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; 28 29 PetscFunctionBegin; 30 ierr = MatGetVariableBlockSizes(pc->pmat,&nblocks,&bsizes);CHKERRQ(ierr); 31 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 32 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 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 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 97 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 98 PetscFunctionReturn(0); 99 } 100 101 102 103 /* -------------------------------------------------------------------------- */ 104 static PetscErrorCode PCSetUp_VPBJacobi(PC pc) 105 { 106 PC_VPBJacobi *jac = (PC_VPBJacobi*)pc->data; 107 PetscErrorCode ierr; 108 Mat A = pc->pmat; 109 MatFactorError err; 110 PetscInt i,nsize = 0,nlocal; 111 PetscInt nblocks; 112 const PetscInt *bsizes; 113 114 PetscFunctionBegin; 115 ierr = MatGetVariableBlockSizes(pc->pmat,&nblocks,&bsizes);CHKERRQ(ierr); 116 ierr = MatGetLocalSize(pc->pmat,&nlocal,NULL);CHKERRQ(ierr); 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 ierr = PetscMalloc1(nsize,&jac->diag);CHKERRQ(ierr); 121 } 122 ierr = MatInvertVariableBlockDiagonal(A,nblocks,bsizes,jac->diag);CHKERRQ(ierr); 123 ierr = MatFactorGetError(A,&err);CHKERRQ(ierr); 124 if (err) pc->failedreason = (PCFailedReason)err; 125 pc->ops->apply = PCApply_VPBJacobi; 126 PetscFunctionReturn(0); 127 } 128 /* -------------------------------------------------------------------------- */ 129 static PetscErrorCode PCDestroy_VPBJacobi(PC pc) 130 { 131 PC_VPBJacobi *jac = (PC_VPBJacobi*)pc->data; 132 PetscErrorCode ierr; 133 134 PetscFunctionBegin; 135 /* 136 Free the private data structure that was hanging off the PC 137 */ 138 ierr = PetscFree(jac->diag);CHKERRQ(ierr); 139 ierr = PetscFree(pc->data);CHKERRQ(ierr); 140 PetscFunctionReturn(0); 141 } 142 143 /* -------------------------------------------------------------------------- */ 144 /*MC 145 PCVPBJACOBI - Variable size point block Jacobi preconditioner 146 147 148 Notes: 149 See PCJACOBI for point Jacobi preconditioning, PCPBJACOBI for fixed point block size, and PCBJACOBI for large size blocks 150 151 This works for AIJ matrices 152 153 Uses dense LU factorization with partial pivoting to invert the blocks; if a zero pivot 154 is detected a PETSc error is generated. 155 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. 162 163 Perhaps should provide an option that allows generation of a valid preconditioner 164 even if a block is singular as the PCJACOBI does. 165 166 Level: beginner 167 168 .seealso: MatSetVariableBlockSizes(), PCCreate(), PCSetType(), PCType (for list of available types), PC, PCJACOBI 169 170 M*/ 171 172 PETSC_EXTERN PetscErrorCode PCCreate_VPBJacobi(PC pc) 173 { 174 PC_VPBJacobi *jac; 175 PetscErrorCode ierr; 176 177 PetscFunctionBegin; 178 /* 179 Creates the private data structure for this preconditioner and 180 attach it to the PC object. 181 */ 182 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr); 183 pc->data = (void*)jac; 184 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; 190 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 PetscFunctionReturn(0); 207 } 208 209 210