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