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