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,*xx,*yy; 49 50 PetscFunctionBegin; 51 ierr = VecGetArray(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 = VecRestoreArray(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,*xx,*yy; 73 74 PetscFunctionBegin; 75 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 76 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 77 for (i=0; i<m; i++) { 78 x0 = xx[3*i]; x1 = xx[3*i+1]; x2 = xx[3*i+2]; 79 80 yy[3*i] = diag[0]*x0 + diag[3]*x1 + diag[6]*x2; 81 yy[3*i+1] = diag[1]*x0 + diag[4]*x1 + diag[7]*x2; 82 yy[3*i+2] = diag[2]*x0 + diag[5]*x1 + diag[8]*x2; 83 diag += 9; 84 } 85 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 86 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 87 ierr = PetscLogFlops(15.0*m);CHKERRQ(ierr); 88 PetscFunctionReturn(0); 89 } 90 #undef __FUNCT__ 91 #define __FUNCT__ "PCApply_PBJacobi_4" 92 static PetscErrorCode PCApply_PBJacobi_4(PC pc,Vec x,Vec y) 93 { 94 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 95 PetscErrorCode ierr; 96 PetscInt i,m = jac->mbs; 97 const MatScalar *diag = jac->diag; 98 PetscScalar x0,x1,x2,x3,*xx,*yy; 99 100 PetscFunctionBegin; 101 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 102 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 103 for (i=0; i<m; i++) { 104 x0 = xx[4*i]; x1 = xx[4*i+1]; x2 = xx[4*i+2]; x3 = xx[4*i+3]; 105 106 yy[4*i] = diag[0]*x0 + diag[4]*x1 + diag[8]*x2 + diag[12]*x3; 107 yy[4*i+1] = diag[1]*x0 + diag[5]*x1 + diag[9]*x2 + diag[13]*x3; 108 yy[4*i+2] = diag[2]*x0 + diag[6]*x1 + diag[10]*x2 + diag[14]*x3; 109 yy[4*i+3] = diag[3]*x0 + diag[7]*x1 + diag[11]*x2 + diag[15]*x3; 110 diag += 16; 111 } 112 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 113 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 114 ierr = PetscLogFlops(28.0*m);CHKERRQ(ierr); 115 PetscFunctionReturn(0); 116 } 117 #undef __FUNCT__ 118 #define __FUNCT__ "PCApply_PBJacobi_5" 119 static PetscErrorCode PCApply_PBJacobi_5(PC pc,Vec x,Vec y) 120 { 121 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 122 PetscErrorCode ierr; 123 PetscInt i,m = jac->mbs; 124 const MatScalar *diag = jac->diag; 125 PetscScalar x0,x1,x2,x3,x4,*xx,*yy; 126 127 PetscFunctionBegin; 128 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 129 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 130 for (i=0; i<m; i++) { 131 x0 = xx[5*i]; x1 = xx[5*i+1]; x2 = xx[5*i+2]; x3 = xx[5*i+3]; x4 = xx[5*i+4]; 132 133 yy[5*i] = diag[0]*x0 + diag[5]*x1 + diag[10]*x2 + diag[15]*x3 + diag[20]*x4; 134 yy[5*i+1] = diag[1]*x0 + diag[6]*x1 + diag[11]*x2 + diag[16]*x3 + diag[21]*x4; 135 yy[5*i+2] = diag[2]*x0 + diag[7]*x1 + diag[12]*x2 + diag[17]*x3 + diag[22]*x4; 136 yy[5*i+3] = diag[3]*x0 + diag[8]*x1 + diag[13]*x2 + diag[18]*x3 + diag[23]*x4; 137 yy[5*i+4] = diag[4]*x0 + diag[9]*x1 + diag[14]*x2 + diag[19]*x3 + diag[24]*x4; 138 diag += 25; 139 } 140 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 141 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 142 ierr = PetscLogFlops(45.0*m);CHKERRQ(ierr); 143 PetscFunctionReturn(0); 144 } 145 #undef __FUNCT__ 146 #define __FUNCT__ "PCApply_PBJacobi_6" 147 static PetscErrorCode PCApply_PBJacobi_6(PC pc,Vec x,Vec y) 148 { 149 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 150 PetscErrorCode ierr; 151 PetscInt i,m = jac->mbs; 152 const MatScalar *diag = jac->diag; 153 PetscScalar x0,x1,x2,x3,x4,x5,*xx,*yy; 154 155 PetscFunctionBegin; 156 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 157 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 158 for (i=0; i<m; i++) { 159 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]; 160 161 yy[6*i] = diag[0]*x0 + diag[6]*x1 + diag[12]*x2 + diag[18]*x3 + diag[24]*x4 + diag[30]*x5; 162 yy[6*i+1] = diag[1]*x0 + diag[7]*x1 + diag[13]*x2 + diag[19]*x3 + diag[25]*x4 + diag[31]*x5; 163 yy[6*i+2] = diag[2]*x0 + diag[8]*x1 + diag[14]*x2 + diag[20]*x3 + diag[26]*x4 + diag[32]*x5; 164 yy[6*i+3] = diag[3]*x0 + diag[9]*x1 + diag[15]*x2 + diag[21]*x3 + diag[27]*x4 + diag[33]*x5; 165 yy[6*i+4] = diag[4]*x0 + diag[10]*x1 + diag[16]*x2 + diag[22]*x3 + diag[28]*x4 + diag[34]*x5; 166 yy[6*i+5] = diag[5]*x0 + diag[11]*x1 + diag[17]*x2 + diag[23]*x3 + diag[29]*x4 + diag[35]*x5; 167 diag += 36; 168 } 169 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 170 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 171 ierr = PetscLogFlops(66.0*m);CHKERRQ(ierr); 172 PetscFunctionReturn(0); 173 } 174 #undef __FUNCT__ 175 #define __FUNCT__ "PCApply_PBJacobi_7" 176 static PetscErrorCode PCApply_PBJacobi_7(PC pc,Vec x,Vec y) 177 { 178 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 179 PetscErrorCode ierr; 180 PetscInt i,m = jac->mbs; 181 const MatScalar *diag = jac->diag; 182 PetscScalar x0,x1,x2,x3,x4,x5,x6,*xx,*yy; 183 184 PetscFunctionBegin; 185 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 186 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 187 for (i=0; i<m; i++) { 188 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]; 189 190 yy[7*i] = diag[0]*x0 + diag[7]*x1 + diag[14]*x2 + diag[21]*x3 + diag[28]*x4 + diag[35]*x5 + diag[42]*x6; 191 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; 192 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; 193 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; 194 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; 195 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; 196 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; 197 diag += 49; 198 } 199 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 200 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 201 ierr = PetscLogFlops(80.0*m);CHKERRQ(ierr); 202 PetscFunctionReturn(0); 203 } 204 /* -------------------------------------------------------------------------- */ 205 #undef __FUNCT__ 206 #define __FUNCT__ "PCSetUp_PBJacobi" 207 static PetscErrorCode PCSetUp_PBJacobi(PC pc) 208 { 209 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 210 PetscErrorCode ierr; 211 Mat A = pc->pmat; 212 213 PetscFunctionBegin; 214 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Supported only for square matrices and square storage"); 215 216 ierr = MatInvertBlockDiagonal(A,&jac->diag);CHKERRQ(ierr); 217 ierr = MatGetBlockSize(A,&jac->bs);CHKERRQ(ierr); 218 jac->mbs = A->rmap->n/jac->bs; 219 switch (jac->bs) { 220 case 1: 221 pc->ops->apply = PCApply_PBJacobi_1; 222 break; 223 case 2: 224 pc->ops->apply = PCApply_PBJacobi_2; 225 break; 226 case 3: 227 pc->ops->apply = PCApply_PBJacobi_3; 228 break; 229 case 4: 230 pc->ops->apply = PCApply_PBJacobi_4; 231 break; 232 case 5: 233 pc->ops->apply = PCApply_PBJacobi_5; 234 break; 235 case 6: 236 pc->ops->apply = PCApply_PBJacobi_6; 237 break; 238 case 7: 239 pc->ops->apply = PCApply_PBJacobi_7; 240 break; 241 default: 242 SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"not supported for block size %D",jac->bs); 243 } 244 PetscFunctionReturn(0); 245 } 246 /* -------------------------------------------------------------------------- */ 247 #undef __FUNCT__ 248 #define __FUNCT__ "PCDestroy_PBJacobi" 249 static PetscErrorCode PCDestroy_PBJacobi(PC pc) 250 { 251 PetscErrorCode ierr; 252 253 PetscFunctionBegin; 254 /* 255 Free the private data structure that was hanging off the PC 256 */ 257 ierr = PetscFree(pc->data);CHKERRQ(ierr); 258 PetscFunctionReturn(0); 259 } 260 261 #undef __FUNCT__ 262 #define __FUNCT__ "PCView_PBJacobi" 263 static PetscErrorCode PCView_PBJacobi(PC pc,PetscViewer viewer) 264 { 265 PetscErrorCode ierr; 266 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 267 PetscBool iascii; 268 269 PetscFunctionBegin; 270 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 271 if (iascii) { 272 ierr = PetscViewerASCIIPrintf(viewer," point-block Jacobi: block size %D\n",jac->bs);CHKERRQ(ierr); 273 } 274 PetscFunctionReturn(0); 275 } 276 277 /* -------------------------------------------------------------------------- */ 278 /*MC 279 PCPBJACOBI - Point block Jacobi 280 281 Level: beginner 282 283 Concepts: point block Jacobi 284 285 286 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC 287 288 M*/ 289 290 #undef __FUNCT__ 291 #define __FUNCT__ "PCCreate_PBJacobi" 292 PETSC_EXTERN PetscErrorCode PCCreate_PBJacobi(PC pc) 293 { 294 PC_PBJacobi *jac; 295 PetscErrorCode ierr; 296 297 PetscFunctionBegin; 298 /* 299 Creates the private data structure for this preconditioner and 300 attach it to the PC object. 301 */ 302 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr); 303 pc->data = (void*)jac; 304 305 /* 306 Initialize the pointers to vectors to ZERO; these will be used to store 307 diagonal entries of the matrix for fast preconditioner application. 308 */ 309 jac->diag = 0; 310 311 /* 312 Set the pointers for the functions that are provided above. 313 Now when the user-level routines (such as PCApply(), PCDestroy(), etc.) 314 are called, they will automatically call these functions. Note we 315 choose not to provide a couple of these functions since they are 316 not needed. 317 */ 318 pc->ops->apply = 0; /*set depending on the block size */ 319 pc->ops->applytranspose = 0; 320 pc->ops->setup = PCSetUp_PBJacobi; 321 pc->ops->destroy = PCDestroy_PBJacobi; 322 pc->ops->setfromoptions = 0; 323 pc->ops->view = PCView_PBJacobi; 324 pc->ops->applyrichardson = 0; 325 pc->ops->applysymmetricleft = 0; 326 pc->ops->applysymmetricright = 0; 327 PetscFunctionReturn(0); 328 } 329 330 331