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 219 PetscFunctionBegin; 220 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Supported only for square matrices and square storage"); 221 222 ierr = MatInvertBlockDiagonal(A,&jac->diag);CHKERRQ(ierr); 223 ierr = MatGetBlockSize(A,&jac->bs);CHKERRQ(ierr); 224 jac->mbs = A->rmap->n/jac->bs; 225 switch (jac->bs) { 226 case 1: 227 pc->ops->apply = PCApply_PBJacobi_1; 228 break; 229 case 2: 230 pc->ops->apply = PCApply_PBJacobi_2; 231 break; 232 case 3: 233 pc->ops->apply = PCApply_PBJacobi_3; 234 break; 235 case 4: 236 pc->ops->apply = PCApply_PBJacobi_4; 237 break; 238 case 5: 239 pc->ops->apply = PCApply_PBJacobi_5; 240 break; 241 case 6: 242 pc->ops->apply = PCApply_PBJacobi_6; 243 break; 244 case 7: 245 pc->ops->apply = PCApply_PBJacobi_7; 246 break; 247 default: 248 SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"not supported for block size %D",jac->bs); 249 } 250 PetscFunctionReturn(0); 251 } 252 /* -------------------------------------------------------------------------- */ 253 #undef __FUNCT__ 254 #define __FUNCT__ "PCDestroy_PBJacobi" 255 static PetscErrorCode PCDestroy_PBJacobi(PC pc) 256 { 257 PetscErrorCode ierr; 258 259 PetscFunctionBegin; 260 /* 261 Free the private data structure that was hanging off the PC 262 */ 263 ierr = PetscFree(pc->data);CHKERRQ(ierr); 264 PetscFunctionReturn(0); 265 } 266 267 #undef __FUNCT__ 268 #define __FUNCT__ "PCView_PBJacobi" 269 static PetscErrorCode PCView_PBJacobi(PC pc,PetscViewer viewer) 270 { 271 PetscErrorCode ierr; 272 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 273 PetscBool iascii; 274 275 PetscFunctionBegin; 276 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 277 if (iascii) { 278 ierr = PetscViewerASCIIPrintf(viewer," point-block Jacobi: block size %D\n",jac->bs);CHKERRQ(ierr); 279 } 280 PetscFunctionReturn(0); 281 } 282 283 /* -------------------------------------------------------------------------- */ 284 /*MC 285 PCPBJACOBI - Point block Jacobi preconditioner 286 287 288 Notes: See PCJACOBI for point Jacobi preconditioning 289 290 This works for AIJ and BAIJ matrices and uses the blocksize provided to the matrix 291 292 Uses dense LU factorization with partial pivoting to invert the blocks; if a zero pivot 293 is detected a PETSc error is generated. 294 295 Developer Notes: This should support the PCSetErrorIfFailure() flag set to PETSC_TRUE to allow 296 the factorization to continue even after a zero pivot is found resulting in a Nan and hence 297 terminating KSP with a KSP_DIVERGED_NANORIF allowing 298 a nonlinear solver/ODE integrator to recover without stopping the program as currently happens. 299 300 Developer Note: Perhaps should provide an option that allows generation of a valid preconditioner 301 even if a block is singular as the PCJACOBI does. 302 303 Level: beginner 304 305 Concepts: point block Jacobi 306 307 308 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCJACOBI 309 310 M*/ 311 312 #undef __FUNCT__ 313 #define __FUNCT__ "PCCreate_PBJacobi" 314 PETSC_EXTERN PetscErrorCode PCCreate_PBJacobi(PC pc) 315 { 316 PC_PBJacobi *jac; 317 PetscErrorCode ierr; 318 319 PetscFunctionBegin; 320 /* 321 Creates the private data structure for this preconditioner and 322 attach it to the PC object. 323 */ 324 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr); 325 pc->data = (void*)jac; 326 327 /* 328 Initialize the pointers to vectors to ZERO; these will be used to store 329 diagonal entries of the matrix for fast preconditioner application. 330 */ 331 jac->diag = 0; 332 333 /* 334 Set the pointers for the functions that are provided above. 335 Now when the user-level routines (such as PCApply(), PCDestroy(), etc.) 336 are called, they will automatically call these functions. Note we 337 choose not to provide a couple of these functions since they are 338 not needed. 339 */ 340 pc->ops->apply = 0; /*set depending on the block size */ 341 pc->ops->applytranspose = 0; 342 pc->ops->setup = PCSetUp_PBJacobi; 343 pc->ops->destroy = PCDestroy_PBJacobi; 344 pc->ops->setfromoptions = 0; 345 pc->ops->view = PCView_PBJacobi; 346 pc->ops->applyrichardson = 0; 347 pc->ops->applysymmetricleft = 0; 348 pc->ops->applysymmetricright = 0; 349 PetscFunctionReturn(0); 350 } 351 352 353