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 /* -------------------------------------------------------------------------- */ 145 #undef __FUNCT__ 146 #define __FUNCT__ "PCSetUp_PBJacobi" 147 static PetscErrorCode PCSetUp_PBJacobi(PC pc) 148 { 149 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 150 PetscErrorCode ierr; 151 Mat A = pc->pmat; 152 153 PetscFunctionBegin; 154 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Supported only for square matrices and square storage"); 155 156 ierr = MatInvertBlockDiagonal(A,&jac->diag);CHKERRQ(ierr); 157 jac->bs = A->rmap->bs; 158 jac->mbs = A->rmap->n/A->rmap->bs; 159 switch (jac->bs){ 160 case 1: 161 pc->ops->apply = PCApply_PBJacobi_1; 162 break; 163 case 2: 164 pc->ops->apply = PCApply_PBJacobi_2; 165 break; 166 case 3: 167 pc->ops->apply = PCApply_PBJacobi_3; 168 break; 169 case 4: 170 pc->ops->apply = PCApply_PBJacobi_4; 171 break; 172 case 5: 173 pc->ops->apply = PCApply_PBJacobi_5; 174 break; 175 default: 176 SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"not supported for block size %D",jac->bs); 177 } 178 179 PetscFunctionReturn(0); 180 } 181 /* -------------------------------------------------------------------------- */ 182 #undef __FUNCT__ 183 #define __FUNCT__ "PCDestroy_PBJacobi" 184 static PetscErrorCode PCDestroy_PBJacobi(PC pc) 185 { 186 PetscErrorCode ierr; 187 188 PetscFunctionBegin; 189 /* 190 Free the private data structure that was hanging off the PC 191 */ 192 ierr = PetscFree(pc->data);CHKERRQ(ierr); 193 PetscFunctionReturn(0); 194 } 195 /* -------------------------------------------------------------------------- */ 196 /*MC 197 PCPBJACOBI - Point block Jacobi 198 199 Level: beginner 200 201 Concepts: point block Jacobi 202 203 204 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC 205 206 M*/ 207 208 EXTERN_C_BEGIN 209 #undef __FUNCT__ 210 #define __FUNCT__ "PCCreate_PBJacobi" 211 PetscErrorCode PCCreate_PBJacobi(PC pc) 212 { 213 PC_PBJacobi *jac; 214 PetscErrorCode ierr; 215 216 PetscFunctionBegin; 217 218 /* 219 Creates the private data structure for this preconditioner and 220 attach it to the PC object. 221 */ 222 ierr = PetscNewLog(pc,PC_PBJacobi,&jac);CHKERRQ(ierr); 223 pc->data = (void*)jac; 224 225 /* 226 Initialize the pointers to vectors to ZERO; these will be used to store 227 diagonal entries of the matrix for fast preconditioner application. 228 */ 229 jac->diag = 0; 230 231 /* 232 Set the pointers for the functions that are provided above. 233 Now when the user-level routines (such as PCApply(), PCDestroy(), etc.) 234 are called, they will automatically call these functions. Note we 235 choose not to provide a couple of these functions since they are 236 not needed. 237 */ 238 pc->ops->apply = 0; /*set depending on the block size */ 239 pc->ops->applytranspose = 0; 240 pc->ops->setup = PCSetUp_PBJacobi; 241 pc->ops->destroy = PCDestroy_PBJacobi; 242 pc->ops->setfromoptions = 0; 243 pc->ops->view = 0; 244 pc->ops->applyrichardson = 0; 245 pc->ops->applysymmetricleft = 0; 246 pc->ops->applysymmetricright = 0; 247 PetscFunctionReturn(0); 248 } 249 EXTERN_C_END 250 251 252