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/pcimpl.h> /*I "petscpc.h" I*/ 8 9 /* 10 Private context (data structure) for the PBJacobi preconditioner. 11 */ 12 typedef struct { 13 MatScalar *diag; 14 PetscInt bs,mbs; 15 } PC_PBJacobi; 16 17 /* 18 Currently only implemented for baij matrices and directly access baij 19 data structures. 20 */ 21 #include <../src/mat/impls/baij/mpi/mpibaij.h> 22 #include <../src/mat/blockinvert.h> 23 24 #undef __FUNCT__ 25 #define __FUNCT__ "PCApply_PBJacobi_1" 26 static PetscErrorCode PCApply_PBJacobi_1(PC pc,Vec x,Vec y) 27 { 28 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 29 PetscErrorCode ierr; 30 PetscInt i,m = jac->mbs; 31 const MatScalar *diag = jac->diag; 32 const PetscScalar *xx; 33 PetscScalar *yy; 34 35 PetscFunctionBegin; 36 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 37 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 38 for (i=0; i<m; i++) { 39 yy[i] = diag[i]*xx[i]; 40 } 41 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 42 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 43 ierr = PetscLogFlops(2.0*m);CHKERRQ(ierr); 44 PetscFunctionReturn(0); 45 } 46 47 #undef __FUNCT__ 48 #define __FUNCT__ "PCApply_PBJacobi_2" 49 static PetscErrorCode PCApply_PBJacobi_2(PC pc,Vec x,Vec y) 50 { 51 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 52 PetscErrorCode ierr; 53 PetscInt i,m = jac->mbs; 54 const MatScalar *diag = jac->diag; 55 PetscScalar x0,x1,*xx,*yy; 56 57 PetscFunctionBegin; 58 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 59 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 60 for (i=0; i<m; i++) { 61 x0 = xx[2*i]; x1 = xx[2*i+1]; 62 yy[2*i] = diag[0]*x0 + diag[2]*x1; 63 yy[2*i+1] = diag[1]*x0 + diag[3]*x1; 64 diag += 4; 65 } 66 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 67 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 68 ierr = PetscLogFlops(6.0*m);CHKERRQ(ierr); 69 PetscFunctionReturn(0); 70 } 71 #undef __FUNCT__ 72 #define __FUNCT__ "PCApply_PBJacobi_3" 73 static PetscErrorCode PCApply_PBJacobi_3(PC pc,Vec x,Vec y) 74 { 75 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 76 PetscErrorCode ierr; 77 PetscInt i,m = jac->mbs; 78 const MatScalar *diag = jac->diag; 79 PetscScalar x0,x1,x2,*xx,*yy; 80 81 PetscFunctionBegin; 82 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 83 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 84 for (i=0; i<m; i++) { 85 x0 = xx[3*i]; x1 = xx[3*i+1]; x2 = xx[3*i+2]; 86 yy[3*i] = diag[0]*x0 + diag[3]*x1 + diag[6]*x2; 87 yy[3*i+1] = diag[1]*x0 + diag[4]*x1 + diag[7]*x2; 88 yy[3*i+2] = diag[2]*x0 + diag[5]*x1 + diag[8]*x2; 89 diag += 9; 90 } 91 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 92 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 93 ierr = PetscLogFlops(15.0*m);CHKERRQ(ierr); 94 PetscFunctionReturn(0); 95 } 96 #undef __FUNCT__ 97 #define __FUNCT__ "PCApply_PBJacobi_4" 98 static PetscErrorCode PCApply_PBJacobi_4(PC pc,Vec x,Vec y) 99 { 100 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 101 PetscErrorCode ierr; 102 PetscInt i,m = jac->mbs; 103 const MatScalar *diag = jac->diag; 104 PetscScalar x0,x1,x2,x3,*xx,*yy; 105 106 PetscFunctionBegin; 107 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 108 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 109 for (i=0; i<m; i++) { 110 x0 = xx[4*i]; x1 = xx[4*i+1]; x2 = xx[4*i+2]; x3 = xx[4*i+3]; 111 yy[4*i] = diag[0]*x0 + diag[4]*x1 + diag[8]*x2 + diag[12]*x3; 112 yy[4*i+1] = diag[1]*x0 + diag[5]*x1 + diag[9]*x2 + diag[13]*x3; 113 yy[4*i+2] = diag[2]*x0 + diag[6]*x1 + diag[10]*x2 + diag[14]*x3; 114 yy[4*i+3] = diag[3]*x0 + diag[7]*x1 + diag[11]*x2 + diag[15]*x3; 115 diag += 16; 116 } 117 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 118 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 119 ierr = PetscLogFlops(28.0*m);CHKERRQ(ierr); 120 PetscFunctionReturn(0); 121 } 122 #undef __FUNCT__ 123 #define __FUNCT__ "PCApply_PBJacobi_5" 124 static PetscErrorCode PCApply_PBJacobi_5(PC pc,Vec x,Vec y) 125 { 126 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 127 PetscErrorCode ierr; 128 PetscInt i,m = jac->mbs; 129 const MatScalar *diag = jac->diag; 130 PetscScalar x0,x1,x2,x3,x4,*xx,*yy; 131 132 PetscFunctionBegin; 133 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 134 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 135 for (i=0; i<m; i++) { 136 x0 = xx[5*i]; x1 = xx[5*i+1]; x2 = xx[5*i+2]; x3 = xx[5*i+3]; x4 = xx[5*i+4]; 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 = VecRestoreArray(x,&xx);CHKERRQ(ierr); 145 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 146 ierr = PetscLogFlops(45.0*m);CHKERRQ(ierr); 147 PetscFunctionReturn(0); 148 } 149 /* -------------------------------------------------------------------------- */ 150 #undef __FUNCT__ 151 #define __FUNCT__ "PCSetUp_PBJacobi" 152 static PetscErrorCode PCSetUp_PBJacobi(PC pc) 153 { 154 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 155 PetscErrorCode ierr; 156 Mat A = pc->pmat; 157 158 PetscFunctionBegin; 159 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Supported only for square matrices and square storage"); 160 161 ierr = MatInvertBlockDiagonal(A,&jac->diag);CHKERRQ(ierr); 162 jac->bs = A->rmap->bs; 163 jac->mbs = A->rmap->n/A->rmap->bs; 164 switch (jac->bs){ 165 case 1: 166 pc->ops->apply = PCApply_PBJacobi_1; 167 break; 168 case 2: 169 pc->ops->apply = PCApply_PBJacobi_2; 170 break; 171 case 3: 172 pc->ops->apply = PCApply_PBJacobi_3; 173 break; 174 case 4: 175 pc->ops->apply = PCApply_PBJacobi_4; 176 break; 177 case 5: 178 pc->ops->apply = PCApply_PBJacobi_5; 179 break; 180 default: 181 SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"not supported for block size %D",jac->bs); 182 } 183 184 PetscFunctionReturn(0); 185 } 186 /* -------------------------------------------------------------------------- */ 187 #undef __FUNCT__ 188 #define __FUNCT__ "PCDestroy_PBJacobi" 189 static PetscErrorCode PCDestroy_PBJacobi(PC pc) 190 { 191 PetscErrorCode ierr; 192 193 PetscFunctionBegin; 194 /* 195 Free the private data structure that was hanging off the PC 196 */ 197 ierr = PetscFree(pc->data);CHKERRQ(ierr); 198 PetscFunctionReturn(0); 199 } 200 /* -------------------------------------------------------------------------- */ 201 /*MC 202 PCPBJACOBI - Point block Jacobi 203 204 Level: beginner 205 206 Concepts: point block Jacobi 207 208 Notes: Only implemented for the BAIJ matrix formats. 209 210 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC 211 212 M*/ 213 214 EXTERN_C_BEGIN 215 #undef __FUNCT__ 216 #define __FUNCT__ "PCCreate_PBJacobi" 217 PetscErrorCode PCCreate_PBJacobi(PC pc) 218 { 219 PC_PBJacobi *jac; 220 PetscErrorCode ierr; 221 222 PetscFunctionBegin; 223 224 /* 225 Creates the private data structure for this preconditioner and 226 attach it to the PC object. 227 */ 228 ierr = PetscNewLog(pc,PC_PBJacobi,&jac);CHKERRQ(ierr); 229 pc->data = (void*)jac; 230 231 /* 232 Initialize the pointers to vectors to ZERO; these will be used to store 233 diagonal entries of the matrix for fast preconditioner application. 234 */ 235 jac->diag = 0; 236 237 /* 238 Set the pointers for the functions that are provided above. 239 Now when the user-level routines (such as PCApply(), PCDestroy(), etc.) 240 are called, they will automatically call these functions. Note we 241 choose not to provide a couple of these functions since they are 242 not needed. 243 */ 244 pc->ops->apply = 0; /*set depending on the block size */ 245 pc->ops->applytranspose = 0; 246 pc->ops->setup = PCSetUp_PBJacobi; 247 pc->ops->destroy = PCDestroy_PBJacobi; 248 pc->ops->setfromoptions = 0; 249 pc->ops->view = 0; 250 pc->ops->applyrichardson = 0; 251 pc->ops->applysymmetricleft = 0; 252 pc->ops->applysymmetricright = 0; 253 PetscFunctionReturn(0); 254 } 255 EXTERN_C_END 256 257 258