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 Concepts: variable point block Jacobi 169 170 .seealso: MatSetVariableBlockSizes(), PCCreate(), PCSetType(), PCType (for list of available types), PC, PCJACOBI 171 172 M*/ 173 174 PETSC_EXTERN PetscErrorCode PCCreate_VPBJacobi(PC pc) 175 { 176 PC_VPBJacobi *jac; 177 PetscErrorCode ierr; 178 179 PetscFunctionBegin; 180 /* 181 Creates the private data structure for this preconditioner and 182 attach it to the PC object. 183 */ 184 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr); 185 pc->data = (void*)jac; 186 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; 192 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 PetscFunctionReturn(0); 209 } 210 211 212