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 PetscInt 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 PetscInt 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 ierr = PetscLogFlops(6*m);CHKERRQ(ierr); 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 PetscInt 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 ierr = PetscLogFlops(15*m);CHKERRQ(ierr); 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 PetscInt 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 ierr = PetscLogFlops(28*m);CHKERRQ(ierr); 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 PetscInt 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 ierr = PetscLogFlops(45*m);CHKERRQ(ierr); 119 PetscFunctionReturn(0); 120 } 121 /* -------------------------------------------------------------------------- */ 122 #undef __FUNCT__ 123 #define __FUNCT__ "PCSetUp_PBJacobi" 124 static PetscErrorCode PCSetUp_PBJacobi(PC pc) 125 { 126 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 127 PetscErrorCode ierr; 128 PetscMPIInt size; 129 PetscTruth seqbaij,mpibaij,baij; 130 Mat A = pc->pmat; 131 Mat_SeqBAIJ *a; 132 133 PetscFunctionBegin; 134 ierr = PetscTypeCompare((PetscObject)pc->pmat,MATSEQBAIJ,&seqbaij);CHKERRQ(ierr); 135 ierr = PetscTypeCompare((PetscObject)pc->pmat,MATMPIBAIJ,&mpibaij);CHKERRQ(ierr); 136 ierr = PetscTypeCompare((PetscObject)pc->pmat,MATBAIJ,&baij);CHKERRQ(ierr); 137 if (!seqbaij && !mpibaij && !baij) { 138 SETERRQ(PETSC_ERR_SUP,"Currently only supports BAIJ matrices"); 139 } 140 ierr = MPI_Comm_size(pc->comm,&size);CHKERRQ(ierr); 141 if (mpibaij || (baij && (size > 1))) A = ((Mat_MPIBAIJ*)A->data)->A; 142 if (A->m != A->n) SETERRQ(PETSC_ERR_SUP,"Supported only for square matrices and square storage"); 143 144 ierr = MatSeqBAIJInvertBlockDiagonal(A);CHKERRQ(ierr); 145 a = (Mat_SeqBAIJ*)A->data; 146 jac->diag = a->idiag; 147 jac->bs = A->bs; 148 jac->mbs = a->mbs; 149 switch (jac->bs){ 150 case 2: 151 pc->ops->apply = PCApply_PBJacobi_2; 152 break; 153 case 3: 154 pc->ops->apply = PCApply_PBJacobi_3; 155 break; 156 case 4: 157 pc->ops->apply = PCApply_PBJacobi_4; 158 break; 159 case 5: 160 pc->ops->apply = PCApply_PBJacobi_5; 161 break; 162 default: 163 SETERRQ1(PETSC_ERR_SUP,"not supported for block size %D",jac->bs); 164 } 165 166 PetscFunctionReturn(0); 167 } 168 /* -------------------------------------------------------------------------- */ 169 #undef __FUNCT__ 170 #define __FUNCT__ "PCDestroy_PBJacobi" 171 static PetscErrorCode PCDestroy_PBJacobi(PC pc) 172 { 173 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 174 PetscErrorCode ierr; 175 176 PetscFunctionBegin; 177 /* 178 Free the private data structure that was hanging off the PC 179 */ 180 ierr = PetscFree(jac);CHKERRQ(ierr); 181 PetscFunctionReturn(0); 182 } 183 /* -------------------------------------------------------------------------- */ 184 /*MC 185 PCPBJACOBI - Point block Jacobi 186 187 Level: beginner 188 189 Concepts: point block Jacobi 190 191 Notes: Only implemented for the BAIJ matrix formats. 192 193 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC 194 195 M*/ 196 197 EXTERN_C_BEGIN 198 #undef __FUNCT__ 199 #define __FUNCT__ "PCCreate_PBJacobi" 200 PetscErrorCode PCCreate_PBJacobi(PC pc) 201 { 202 PC_PBJacobi *jac; 203 PetscErrorCode ierr; 204 205 PetscFunctionBegin; 206 207 /* 208 Creates the private data structure for this preconditioner and 209 attach it to the PC object. 210 */ 211 ierr = PetscNew(PC_PBJacobi,&jac);CHKERRQ(ierr); 212 pc->data = (void*)jac; 213 214 /* 215 Logs the memory usage; this is not needed but allows PETSc to 216 monitor how much memory is being used for various purposes. 217 */ 218 ierr = PetscLogObjectMemory(pc,sizeof(PC_PBJacobi));CHKERRQ(ierr); 219 220 /* 221 Initialize the pointers to vectors to ZERO; these will be used to store 222 diagonal entries of the matrix for fast preconditioner application. 223 */ 224 jac->diag = 0; 225 226 /* 227 Set the pointers for the functions that are provided above. 228 Now when the user-level routines (such as PCApply(), PCDestroy(), etc.) 229 are called, they will automatically call these functions. Note we 230 choose not to provide a couple of these functions since they are 231 not needed. 232 */ 233 pc->ops->apply = 0; /*set depending on the block size */ 234 pc->ops->applytranspose = 0; 235 pc->ops->setup = PCSetUp_PBJacobi; 236 pc->ops->destroy = PCDestroy_PBJacobi; 237 pc->ops->setfromoptions = 0; 238 pc->ops->view = 0; 239 pc->ops->applyrichardson = 0; 240 pc->ops->applysymmetricleft = 0; 241 pc->ops->applysymmetricright = 0; 242 PetscFunctionReturn(0); 243 } 244 EXTERN_C_END 245 246 247