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