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