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 286 287 Level: beginner 288 289 Concepts: point block Jacobi 290 291 292 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC 293 294 M*/ 295 296 #undef __FUNCT__ 297 #define __FUNCT__ "PCCreate_PBJacobi" 298 PETSC_EXTERN PetscErrorCode PCCreate_PBJacobi(PC pc) 299 { 300 PC_PBJacobi *jac; 301 PetscErrorCode ierr; 302 303 PetscFunctionBegin; 304 /* 305 Creates the private data structure for this preconditioner and 306 attach it to the PC object. 307 */ 308 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr); 309 pc->data = (void*)jac; 310 311 /* 312 Initialize the pointers to vectors to ZERO; these will be used to store 313 diagonal entries of the matrix for fast preconditioner application. 314 */ 315 jac->diag = 0; 316 317 /* 318 Set the pointers for the functions that are provided above. 319 Now when the user-level routines (such as PCApply(), PCDestroy(), etc.) 320 are called, they will automatically call these functions. Note we 321 choose not to provide a couple of these functions since they are 322 not needed. 323 */ 324 pc->ops->apply = 0; /*set depending on the block size */ 325 pc->ops->applytranspose = 0; 326 pc->ops->setup = PCSetUp_PBJacobi; 327 pc->ops->destroy = PCDestroy_PBJacobi; 328 pc->ops->setfromoptions = 0; 329 pc->ops->view = PCView_PBJacobi; 330 pc->ops->applyrichardson = 0; 331 pc->ops->applysymmetricleft = 0; 332 pc->ops->applysymmetricright = 0; 333 PetscFunctionReturn(0); 334 } 335 336 337