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 <private/matimpl.h> 8 #include <private/pcimpl.h> /*I "petscpc.h" I*/ 9 10 /* 11 Private context (data structure) for the PBJacobi preconditioner. 12 */ 13 typedef struct { 14 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 /* -------------------------------------------------------------------------- */ 173 #undef __FUNCT__ 174 #define __FUNCT__ "PCSetUp_PBJacobi" 175 static PetscErrorCode PCSetUp_PBJacobi(PC pc) 176 { 177 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 178 PetscErrorCode ierr; 179 Mat A = pc->pmat; 180 181 PetscFunctionBegin; 182 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Supported only for square matrices and square storage"); 183 184 ierr = MatInvertBlockDiagonal(A,&jac->diag);CHKERRQ(ierr); 185 jac->bs = A->rmap->bs; 186 jac->mbs = A->rmap->n/A->rmap->bs; 187 switch (jac->bs){ 188 case 1: 189 pc->ops->apply = PCApply_PBJacobi_1; 190 break; 191 case 2: 192 pc->ops->apply = PCApply_PBJacobi_2; 193 break; 194 case 3: 195 pc->ops->apply = PCApply_PBJacobi_3; 196 break; 197 case 4: 198 pc->ops->apply = PCApply_PBJacobi_4; 199 break; 200 case 5: 201 pc->ops->apply = PCApply_PBJacobi_5; 202 break; 203 case 6: 204 pc->ops->apply = PCApply_PBJacobi_6; 205 break; 206 default: 207 SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"not supported for block size %D",jac->bs); 208 } 209 210 PetscFunctionReturn(0); 211 } 212 /* -------------------------------------------------------------------------- */ 213 #undef __FUNCT__ 214 #define __FUNCT__ "PCDestroy_PBJacobi" 215 static PetscErrorCode PCDestroy_PBJacobi(PC pc) 216 { 217 PetscErrorCode ierr; 218 219 PetscFunctionBegin; 220 /* 221 Free the private data structure that was hanging off the PC 222 */ 223 ierr = PetscFree(pc->data);CHKERRQ(ierr); 224 PetscFunctionReturn(0); 225 } 226 /* -------------------------------------------------------------------------- */ 227 /*MC 228 PCPBJACOBI - Point block Jacobi 229 230 Level: beginner 231 232 Concepts: point block Jacobi 233 234 235 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC 236 237 M*/ 238 239 EXTERN_C_BEGIN 240 #undef __FUNCT__ 241 #define __FUNCT__ "PCCreate_PBJacobi" 242 PetscErrorCode PCCreate_PBJacobi(PC pc) 243 { 244 PC_PBJacobi *jac; 245 PetscErrorCode ierr; 246 247 PetscFunctionBegin; 248 249 /* 250 Creates the private data structure for this preconditioner and 251 attach it to the PC object. 252 */ 253 ierr = PetscNewLog(pc,PC_PBJacobi,&jac);CHKERRQ(ierr); 254 pc->data = (void*)jac; 255 256 /* 257 Initialize the pointers to vectors to ZERO; these will be used to store 258 diagonal entries of the matrix for fast preconditioner application. 259 */ 260 jac->diag = 0; 261 262 /* 263 Set the pointers for the functions that are provided above. 264 Now when the user-level routines (such as PCApply(), PCDestroy(), etc.) 265 are called, they will automatically call these functions. Note we 266 choose not to provide a couple of these functions since they are 267 not needed. 268 */ 269 pc->ops->apply = 0; /*set depending on the block size */ 270 pc->ops->applytranspose = 0; 271 pc->ops->setup = PCSetUp_PBJacobi; 272 pc->ops->destroy = PCDestroy_PBJacobi; 273 pc->ops->setfromoptions = 0; 274 pc->ops->view = 0; 275 pc->ops->applyrichardson = 0; 276 pc->ops->applysymmetricleft = 0; 277 pc->ops->applysymmetricright = 0; 278 PetscFunctionReturn(0); 279 } 280 EXTERN_C_END 281 282 283