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++) { 34 yy[i] = diag[i]*xx[i]; 35 } 36 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 37 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 38 ierr = PetscLogFlops(2.0*m);CHKERRQ(ierr); 39 PetscFunctionReturn(0); 40 } 41 42 #undef __FUNCT__ 43 #define __FUNCT__ "PCApply_PBJacobi_2" 44 static PetscErrorCode PCApply_PBJacobi_2(PC pc,Vec x,Vec y) 45 { 46 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 47 PetscErrorCode ierr; 48 PetscInt i,m = jac->mbs; 49 const MatScalar *diag = jac->diag; 50 PetscScalar x0,x1,*xx,*yy; 51 52 PetscFunctionBegin; 53 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 54 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 55 for (i=0; i<m; i++) { 56 x0 = xx[2*i]; x1 = xx[2*i+1]; 57 yy[2*i] = diag[0]*x0 + diag[2]*x1; 58 yy[2*i+1] = diag[1]*x0 + diag[3]*x1; 59 diag += 4; 60 } 61 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 62 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 63 ierr = PetscLogFlops(6.0*m);CHKERRQ(ierr); 64 PetscFunctionReturn(0); 65 } 66 #undef __FUNCT__ 67 #define __FUNCT__ "PCApply_PBJacobi_3" 68 static PetscErrorCode PCApply_PBJacobi_3(PC pc,Vec x,Vec y) 69 { 70 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 71 PetscErrorCode ierr; 72 PetscInt i,m = jac->mbs; 73 const MatScalar *diag = jac->diag; 74 PetscScalar x0,x1,x2,*xx,*yy; 75 76 PetscFunctionBegin; 77 ierr = VecGetArray(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 yy[3*i] = diag[0]*x0 + diag[3]*x1 + diag[6]*x2; 82 yy[3*i+1] = diag[1]*x0 + diag[4]*x1 + diag[7]*x2; 83 yy[3*i+2] = diag[2]*x0 + diag[5]*x1 + diag[8]*x2; 84 diag += 9; 85 } 86 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 87 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 88 ierr = PetscLogFlops(15.0*m);CHKERRQ(ierr); 89 PetscFunctionReturn(0); 90 } 91 #undef __FUNCT__ 92 #define __FUNCT__ "PCApply_PBJacobi_4" 93 static PetscErrorCode PCApply_PBJacobi_4(PC pc,Vec x,Vec y) 94 { 95 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 96 PetscErrorCode ierr; 97 PetscInt i,m = jac->mbs; 98 const MatScalar *diag = jac->diag; 99 PetscScalar x0,x1,x2,x3,*xx,*yy; 100 101 PetscFunctionBegin; 102 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 103 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 104 for (i=0; i<m; i++) { 105 x0 = xx[4*i]; x1 = xx[4*i+1]; x2 = xx[4*i+2]; x3 = xx[4*i+3]; 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 yy[5*i] = diag[0]*x0 + diag[5]*x1 + diag[10]*x2 + diag[15]*x3 + diag[20]*x4; 133 yy[5*i+1] = diag[1]*x0 + diag[6]*x1 + diag[11]*x2 + diag[16]*x3 + diag[21]*x4; 134 yy[5*i+2] = diag[2]*x0 + diag[7]*x1 + diag[12]*x2 + diag[17]*x3 + diag[22]*x4; 135 yy[5*i+3] = diag[3]*x0 + diag[8]*x1 + diag[13]*x2 + diag[18]*x3 + diag[23]*x4; 136 yy[5*i+4] = diag[4]*x0 + diag[9]*x1 + diag[14]*x2 + diag[19]*x3 + diag[24]*x4; 137 diag += 25; 138 } 139 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 140 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 141 ierr = PetscLogFlops(45.0*m);CHKERRQ(ierr); 142 PetscFunctionReturn(0); 143 } 144 #undef __FUNCT__ 145 #define __FUNCT__ "PCApply_PBJacobi_6" 146 static PetscErrorCode PCApply_PBJacobi_6(PC pc,Vec x,Vec y) 147 { 148 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 149 PetscErrorCode ierr; 150 PetscInt i,m = jac->mbs; 151 const MatScalar *diag = jac->diag; 152 PetscScalar x0,x1,x2,x3,x4,x5,*xx,*yy; 153 154 PetscFunctionBegin; 155 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 156 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 157 for (i=0; i<m; i++) { 158 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]; 159 yy[6*i] = diag[0]*x0 + diag[6]*x1 + diag[12]*x2 + diag[18]*x3 + diag[24]*x4 + diag[30]*x5; 160 yy[6*i+1] = diag[1]*x0 + diag[7]*x1 + diag[13]*x2 + diag[19]*x3 + diag[25]*x4 + diag[31]*x5; 161 yy[6*i+2] = diag[2]*x0 + diag[8]*x1 + diag[14]*x2 + diag[20]*x3 + diag[26]*x4 + diag[32]*x5; 162 yy[6*i+3] = diag[3]*x0 + diag[9]*x1 + diag[15]*x2 + diag[21]*x3 + diag[27]*x4 + diag[33]*x5; 163 yy[6*i+4] = diag[4]*x0 + diag[10]*x1 + diag[16]*x2 + diag[22]*x3 + diag[28]*x4 + diag[34]*x5; 164 yy[6*i+5] = diag[5]*x0 + diag[11]*x1 + diag[17]*x2 + diag[23]*x3 + diag[29]*x4 + diag[35]*x5; 165 diag += 36; 166 } 167 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 168 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 169 ierr = PetscLogFlops(66.0*m);CHKERRQ(ierr); 170 PetscFunctionReturn(0); 171 } 172 #undef __FUNCT__ 173 #define __FUNCT__ "PCApply_PBJacobi_7" 174 static PetscErrorCode PCApply_PBJacobi_7(PC pc,Vec x,Vec y) 175 { 176 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 177 PetscErrorCode ierr; 178 PetscInt i,m = jac->mbs; 179 const MatScalar *diag = jac->diag; 180 PetscScalar x0,x1,x2,x3,x4,x5,x6,*xx,*yy; 181 182 PetscFunctionBegin; 183 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 184 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 185 for (i=0; i<m; i++) { 186 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]; 187 yy[7*i] = diag[0]*x0 + diag[7]*x1 + diag[14]*x2 + diag[21]*x3 + diag[28]*x4 + diag[35]*x5 + diag[42]*x6; 188 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; 189 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; 190 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; 191 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; 192 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; 193 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; 194 diag += 49; 195 } 196 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 197 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 198 ierr = PetscLogFlops(80.0*m);CHKERRQ(ierr); 199 PetscFunctionReturn(0); 200 } 201 /* -------------------------------------------------------------------------- */ 202 #undef __FUNCT__ 203 #define __FUNCT__ "PCSetUp_PBJacobi" 204 static PetscErrorCode PCSetUp_PBJacobi(PC pc) 205 { 206 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 207 PetscErrorCode ierr; 208 Mat A = pc->pmat; 209 210 PetscFunctionBegin; 211 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Supported only for square matrices and square storage"); 212 213 ierr = MatInvertBlockDiagonal(A,&jac->diag);CHKERRQ(ierr); 214 jac->bs = A->rmap->bs; 215 jac->mbs = A->rmap->n/A->rmap->bs; 216 switch (jac->bs){ 217 case 1: 218 pc->ops->apply = PCApply_PBJacobi_1; 219 break; 220 case 2: 221 pc->ops->apply = PCApply_PBJacobi_2; 222 break; 223 case 3: 224 pc->ops->apply = PCApply_PBJacobi_3; 225 break; 226 case 4: 227 pc->ops->apply = PCApply_PBJacobi_4; 228 break; 229 case 5: 230 pc->ops->apply = PCApply_PBJacobi_5; 231 break; 232 case 6: 233 pc->ops->apply = PCApply_PBJacobi_6; 234 break; 235 case 7: 236 pc->ops->apply = PCApply_PBJacobi_7; 237 break; 238 default: 239 SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"not supported for block size %D",jac->bs); 240 } 241 242 PetscFunctionReturn(0); 243 } 244 /* -------------------------------------------------------------------------- */ 245 #undef __FUNCT__ 246 #define __FUNCT__ "PCDestroy_PBJacobi" 247 static PetscErrorCode PCDestroy_PBJacobi(PC pc) 248 { 249 PetscErrorCode ierr; 250 251 PetscFunctionBegin; 252 /* 253 Free the private data structure that was hanging off the PC 254 */ 255 ierr = PetscFree(pc->data);CHKERRQ(ierr); 256 PetscFunctionReturn(0); 257 } 258 259 #undef __FUNCT__ 260 #define __FUNCT__ "PCView_PBJacobi" 261 static PetscErrorCode PCView_PBJacobi(PC pc,PetscViewer viewer) 262 { 263 PetscErrorCode ierr; 264 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 265 PetscBool iascii; 266 267 PetscFunctionBegin; 268 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 269 if (iascii) { 270 ierr = PetscViewerASCIIPrintf(viewer," point-block Jacobi: block size %D\n",jac->bs);CHKERRQ(ierr); 271 } else SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for point-block Jacobi",((PetscObject)viewer)->type_name); 272 PetscFunctionReturn(0); 273 } 274 275 /* -------------------------------------------------------------------------- */ 276 /*MC 277 PCPBJACOBI - Point block Jacobi 278 279 Level: beginner 280 281 Concepts: point block Jacobi 282 283 284 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC 285 286 M*/ 287 288 EXTERN_C_BEGIN 289 #undef __FUNCT__ 290 #define __FUNCT__ "PCCreate_PBJacobi" 291 PetscErrorCode PCCreate_PBJacobi(PC pc) 292 { 293 PC_PBJacobi *jac; 294 PetscErrorCode ierr; 295 296 PetscFunctionBegin; 297 298 /* 299 Creates the private data structure for this preconditioner and 300 attach it to the PC object. 301 */ 302 ierr = PetscNewLog(pc,PC_PBJacobi,&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 EXTERN_C_END 330 331 332