1 2 /* 3 Include files needed for the 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 PBJacobi preconditioner. 11 */ 12 typedef struct { 13 const MatScalar *diag; 14 PetscInt bs,mbs; 15 } PC_PBJacobi; 16 17 18 #undef __FUNCT__ 19 #define __FUNCT__ "PCApply_PBJacobi_1" 20 static PetscErrorCode PCApply_PBJacobi_1(PC pc,Vec x,Vec y) 21 { 22 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 23 PetscErrorCode ierr; 24 PetscInt i,m = jac->mbs; 25 const MatScalar *diag = jac->diag; 26 const PetscScalar *xx; 27 PetscScalar *yy; 28 29 PetscFunctionBegin; 30 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 31 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 32 for (i=0; i<m; i++) yy[i] = diag[i]*xx[i]; 33 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 34 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 35 ierr = PetscLogFlops(2.0*m);CHKERRQ(ierr); 36 PetscFunctionReturn(0); 37 } 38 39 #undef __FUNCT__ 40 #define __FUNCT__ "PCApply_PBJacobi_2" 41 static PetscErrorCode PCApply_PBJacobi_2(PC pc,Vec x,Vec y) 42 { 43 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 44 PetscErrorCode ierr; 45 PetscInt i,m = jac->mbs; 46 const MatScalar *diag = jac->diag; 47 PetscScalar x0,x1,*yy; 48 const PetscScalar *xx; 49 50 PetscFunctionBegin; 51 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 52 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 53 for (i=0; i<m; i++) { 54 x0 = xx[2*i]; x1 = xx[2*i+1]; 55 yy[2*i] = diag[0]*x0 + diag[2]*x1; 56 yy[2*i+1] = diag[1]*x0 + diag[3]*x1; 57 diag += 4; 58 } 59 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 60 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 61 ierr = PetscLogFlops(6.0*m);CHKERRQ(ierr); 62 PetscFunctionReturn(0); 63 } 64 #undef __FUNCT__ 65 #define __FUNCT__ "PCApply_PBJacobi_3" 66 static PetscErrorCode PCApply_PBJacobi_3(PC pc,Vec x,Vec y) 67 { 68 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 69 PetscErrorCode ierr; 70 PetscInt i,m = jac->mbs; 71 const MatScalar *diag = jac->diag; 72 PetscScalar x0,x1,x2,*yy; 73 const PetscScalar *xx; 74 75 PetscFunctionBegin; 76 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 77 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 78 for (i=0; i<m; i++) { 79 x0 = xx[3*i]; x1 = xx[3*i+1]; x2 = xx[3*i+2]; 80 81 yy[3*i] = diag[0]*x0 + diag[3]*x1 + diag[6]*x2; 82 yy[3*i+1] = diag[1]*x0 + diag[4]*x1 + diag[7]*x2; 83 yy[3*i+2] = diag[2]*x0 + diag[5]*x1 + diag[8]*x2; 84 diag += 9; 85 } 86 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 87 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 88 ierr = PetscLogFlops(15.0*m);CHKERRQ(ierr); 89 PetscFunctionReturn(0); 90 } 91 #undef __FUNCT__ 92 #define __FUNCT__ "PCApply_PBJacobi_4" 93 static PetscErrorCode PCApply_PBJacobi_4(PC pc,Vec x,Vec y) 94 { 95 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 96 PetscErrorCode ierr; 97 PetscInt i,m = jac->mbs; 98 const MatScalar *diag = jac->diag; 99 PetscScalar x0,x1,x2,x3,*yy; 100 const PetscScalar *xx; 101 102 PetscFunctionBegin; 103 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 104 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 105 for (i=0; i<m; i++) { 106 x0 = xx[4*i]; x1 = xx[4*i+1]; x2 = xx[4*i+2]; x3 = xx[4*i+3]; 107 108 yy[4*i] = diag[0]*x0 + diag[4]*x1 + diag[8]*x2 + diag[12]*x3; 109 yy[4*i+1] = diag[1]*x0 + diag[5]*x1 + diag[9]*x2 + diag[13]*x3; 110 yy[4*i+2] = diag[2]*x0 + diag[6]*x1 + diag[10]*x2 + diag[14]*x3; 111 yy[4*i+3] = diag[3]*x0 + diag[7]*x1 + diag[11]*x2 + diag[15]*x3; 112 diag += 16; 113 } 114 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 115 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 116 ierr = PetscLogFlops(28.0*m);CHKERRQ(ierr); 117 PetscFunctionReturn(0); 118 } 119 #undef __FUNCT__ 120 #define __FUNCT__ "PCApply_PBJacobi_5" 121 static PetscErrorCode PCApply_PBJacobi_5(PC pc,Vec x,Vec y) 122 { 123 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 124 PetscErrorCode ierr; 125 PetscInt i,m = jac->mbs; 126 const MatScalar *diag = jac->diag; 127 PetscScalar x0,x1,x2,x3,x4,*yy; 128 const PetscScalar *xx; 129 130 PetscFunctionBegin; 131 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 132 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 133 for (i=0; i<m; i++) { 134 x0 = xx[5*i]; x1 = xx[5*i+1]; x2 = xx[5*i+2]; x3 = xx[5*i+3]; x4 = xx[5*i+4]; 135 136 yy[5*i] = diag[0]*x0 + diag[5]*x1 + diag[10]*x2 + diag[15]*x3 + diag[20]*x4; 137 yy[5*i+1] = diag[1]*x0 + diag[6]*x1 + diag[11]*x2 + diag[16]*x3 + diag[21]*x4; 138 yy[5*i+2] = diag[2]*x0 + diag[7]*x1 + diag[12]*x2 + diag[17]*x3 + diag[22]*x4; 139 yy[5*i+3] = diag[3]*x0 + diag[8]*x1 + diag[13]*x2 + diag[18]*x3 + diag[23]*x4; 140 yy[5*i+4] = diag[4]*x0 + diag[9]*x1 + diag[14]*x2 + diag[19]*x3 + diag[24]*x4; 141 diag += 25; 142 } 143 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 144 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 145 ierr = PetscLogFlops(45.0*m);CHKERRQ(ierr); 146 PetscFunctionReturn(0); 147 } 148 #undef __FUNCT__ 149 #define __FUNCT__ "PCApply_PBJacobi_6" 150 static PetscErrorCode PCApply_PBJacobi_6(PC pc,Vec x,Vec y) 151 { 152 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 153 PetscErrorCode ierr; 154 PetscInt i,m = jac->mbs; 155 const MatScalar *diag = jac->diag; 156 PetscScalar x0,x1,x2,x3,x4,x5,*yy; 157 const PetscScalar *xx; 158 159 PetscFunctionBegin; 160 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 161 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 162 for (i=0; i<m; i++) { 163 x0 = xx[6*i]; x1 = xx[6*i+1]; x2 = xx[6*i+2]; x3 = xx[6*i+3]; x4 = xx[6*i+4]; x5 = xx[6*i+5]; 164 165 yy[6*i] = diag[0]*x0 + diag[6]*x1 + diag[12]*x2 + diag[18]*x3 + diag[24]*x4 + diag[30]*x5; 166 yy[6*i+1] = diag[1]*x0 + diag[7]*x1 + diag[13]*x2 + diag[19]*x3 + diag[25]*x4 + diag[31]*x5; 167 yy[6*i+2] = diag[2]*x0 + diag[8]*x1 + diag[14]*x2 + diag[20]*x3 + diag[26]*x4 + diag[32]*x5; 168 yy[6*i+3] = diag[3]*x0 + diag[9]*x1 + diag[15]*x2 + diag[21]*x3 + diag[27]*x4 + diag[33]*x5; 169 yy[6*i+4] = diag[4]*x0 + diag[10]*x1 + diag[16]*x2 + diag[22]*x3 + diag[28]*x4 + diag[34]*x5; 170 yy[6*i+5] = diag[5]*x0 + diag[11]*x1 + diag[17]*x2 + diag[23]*x3 + diag[29]*x4 + diag[35]*x5; 171 diag += 36; 172 } 173 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 174 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 175 ierr = PetscLogFlops(66.0*m);CHKERRQ(ierr); 176 PetscFunctionReturn(0); 177 } 178 #undef __FUNCT__ 179 #define __FUNCT__ "PCApply_PBJacobi_7" 180 static PetscErrorCode PCApply_PBJacobi_7(PC pc,Vec x,Vec y) 181 { 182 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 183 PetscErrorCode ierr; 184 PetscInt i,m = jac->mbs; 185 const MatScalar *diag = jac->diag; 186 PetscScalar x0,x1,x2,x3,x4,x5,x6,*yy; 187 const PetscScalar *xx; 188 189 PetscFunctionBegin; 190 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 191 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 192 for (i=0; i<m; i++) { 193 x0 = xx[7*i]; x1 = xx[7*i+1]; x2 = xx[7*i+2]; x3 = xx[7*i+3]; x4 = xx[7*i+4]; x5 = xx[7*i+5]; x6 = xx[7*i+6]; 194 195 yy[7*i] = diag[0]*x0 + diag[7]*x1 + diag[14]*x2 + diag[21]*x3 + diag[28]*x4 + diag[35]*x5 + diag[42]*x6; 196 yy[7*i+1] = diag[1]*x0 + diag[8]*x1 + diag[15]*x2 + diag[22]*x3 + diag[29]*x4 + diag[36]*x5 + diag[43]*x6; 197 yy[7*i+2] = diag[2]*x0 + diag[9]*x1 + diag[16]*x2 + diag[23]*x3 + diag[30]*x4 + diag[37]*x5 + diag[44]*x6; 198 yy[7*i+3] = diag[3]*x0 + diag[10]*x1 + diag[17]*x2 + diag[24]*x3 + diag[31]*x4 + diag[38]*x5 + diag[45]*x6; 199 yy[7*i+4] = diag[4]*x0 + diag[11]*x1 + diag[18]*x2 + diag[25]*x3 + diag[32]*x4 + diag[39]*x5 + diag[46]*x6; 200 yy[7*i+5] = diag[5]*x0 + diag[12]*x1 + diag[19]*x2 + diag[26]*x3 + diag[33]*x4 + diag[40]*x5 + diag[47]*x6; 201 yy[7*i+6] = diag[6]*x0 + diag[13]*x1 + diag[20]*x2 + diag[27]*x3 + diag[34]*x4 + diag[41]*x5 + diag[48]*x6; 202 diag += 49; 203 } 204 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 205 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 206 ierr = PetscLogFlops(80.0*m);CHKERRQ(ierr); 207 PetscFunctionReturn(0); 208 } 209 /* -------------------------------------------------------------------------- */ 210 #undef __FUNCT__ 211 #define __FUNCT__ "PCSetUp_PBJacobi" 212 static PetscErrorCode PCSetUp_PBJacobi(PC pc) 213 { 214 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 215 PetscErrorCode ierr; 216 Mat A = pc->pmat; 217 MatFactorError err; 218 PetscInt nlocal; 219 220 PetscFunctionBegin; 221 ierr = MatInvertBlockDiagonal(A,&jac->diag);CHKERRQ(ierr); 222 ierr = MatFactorGetError(A,&err);CHKERRQ(ierr); 223 if (err) pc->failedreason = (PCFailedReason)err; 224 225 ierr = MatGetBlockSize(A,&jac->bs);CHKERRQ(ierr); 226 ierr = MatGetLocalSize(A,&nlocal,NULL);CHKERRQ(ierr); 227 jac->mbs = nlocal/jac->bs; 228 switch (jac->bs) { 229 case 1: 230 pc->ops->apply = PCApply_PBJacobi_1; 231 break; 232 case 2: 233 pc->ops->apply = PCApply_PBJacobi_2; 234 break; 235 case 3: 236 pc->ops->apply = PCApply_PBJacobi_3; 237 break; 238 case 4: 239 pc->ops->apply = PCApply_PBJacobi_4; 240 break; 241 case 5: 242 pc->ops->apply = PCApply_PBJacobi_5; 243 break; 244 case 6: 245 pc->ops->apply = PCApply_PBJacobi_6; 246 break; 247 case 7: 248 pc->ops->apply = PCApply_PBJacobi_7; 249 break; 250 default: 251 SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"not supported for block size %D",jac->bs); 252 } 253 PetscFunctionReturn(0); 254 } 255 /* -------------------------------------------------------------------------- */ 256 #undef __FUNCT__ 257 #define __FUNCT__ "PCDestroy_PBJacobi" 258 static PetscErrorCode PCDestroy_PBJacobi(PC pc) 259 { 260 PetscErrorCode ierr; 261 262 PetscFunctionBegin; 263 /* 264 Free the private data structure that was hanging off the PC 265 */ 266 ierr = PetscFree(pc->data);CHKERRQ(ierr); 267 PetscFunctionReturn(0); 268 } 269 270 #undef __FUNCT__ 271 #define __FUNCT__ "PCView_PBJacobi" 272 static PetscErrorCode PCView_PBJacobi(PC pc,PetscViewer viewer) 273 { 274 PetscErrorCode ierr; 275 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 276 PetscBool iascii; 277 278 PetscFunctionBegin; 279 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 280 if (iascii) { 281 ierr = PetscViewerASCIIPrintf(viewer," point-block Jacobi: block size %D\n",jac->bs);CHKERRQ(ierr); 282 } 283 PetscFunctionReturn(0); 284 } 285 286 /* -------------------------------------------------------------------------- */ 287 /*MC 288 PCPBJACOBI - Point block Jacobi preconditioner 289 290 291 Notes: See PCJACOBI for point Jacobi preconditioning 292 293 This works for AIJ and BAIJ matrices and uses the blocksize provided to the matrix 294 295 Uses dense LU factorization with partial pivoting to invert the blocks; if a zero pivot 296 is detected a PETSc error is generated. 297 298 Developer Notes: This should support the PCSetErrorIfFailure() flag set to PETSC_TRUE to allow 299 the factorization to continue even after a zero pivot is found resulting in a Nan and hence 300 terminating KSP with a KSP_DIVERGED_NANORIF allowing 301 a nonlinear solver/ODE integrator to recover without stopping the program as currently happens. 302 303 Developer Note: Perhaps should provide an option that allows generation of a valid preconditioner 304 even if a block is singular as the PCJACOBI does. 305 306 Level: beginner 307 308 Concepts: point block Jacobi 309 310 311 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCJACOBI 312 313 M*/ 314 315 #undef __FUNCT__ 316 #define __FUNCT__ "PCCreate_PBJacobi" 317 PETSC_EXTERN PetscErrorCode PCCreate_PBJacobi(PC pc) 318 { 319 PC_PBJacobi *jac; 320 PetscErrorCode ierr; 321 322 PetscFunctionBegin; 323 /* 324 Creates the private data structure for this preconditioner and 325 attach it to the PC object. 326 */ 327 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr); 328 pc->data = (void*)jac; 329 330 /* 331 Initialize the pointers to vectors to ZERO; these will be used to store 332 diagonal entries of the matrix for fast preconditioner application. 333 */ 334 jac->diag = 0; 335 336 /* 337 Set the pointers for the functions that are provided above. 338 Now when the user-level routines (such as PCApply(), PCDestroy(), etc.) 339 are called, they will automatically call these functions. Note we 340 choose not to provide a couple of these functions since they are 341 not needed. 342 */ 343 pc->ops->apply = 0; /*set depending on the block size */ 344 pc->ops->applytranspose = 0; 345 pc->ops->setup = PCSetUp_PBJacobi; 346 pc->ops->destroy = PCDestroy_PBJacobi; 347 pc->ops->setfromoptions = 0; 348 pc->ops->view = PCView_PBJacobi; 349 pc->ops->applyrichardson = 0; 350 pc->ops->applysymmetricleft = 0; 351 pc->ops->applysymmetricright = 0; 352 PetscFunctionReturn(0); 353 } 354 355 356