1 #define PETSCKSP_DLL 2 3 /* 4 Include files needed for the PBJacobi preconditioner: 5 pcimpl.h - private include file intended for use by all preconditioners 6 */ 7 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 PetscScalar *diag; 15 PetscInt bs,mbs; 16 } PC_PBJacobi; 17 18 /* 19 Currently only implemented for baij matrices and directly access baij 20 data structures. 21 */ 22 #include "src/mat/impls/baij/mpi/mpibaij.h" 23 #include "src/inline/ilu.h" 24 25 #undef __FUNCT__ 26 #define __FUNCT__ "PCApply_PBJacobi_2" 27 static PetscErrorCode PCApply_PBJacobi_2(PC pc,Vec x,Vec y) 28 { 29 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 30 PetscErrorCode ierr; 31 PetscInt i,m = jac->mbs; 32 PetscScalar *diag = jac->diag,x0,x1,*xx,*yy; 33 34 PetscFunctionBegin; 35 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 36 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 37 for (i=0; i<m; i++) { 38 x0 = xx[2*i]; x1 = xx[2*i+1]; 39 yy[2*i] = diag[0]*x0 + diag[2]*x1; 40 yy[2*i+1] = diag[1]*x0 + diag[3]*x1; 41 diag += 4; 42 } 43 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 44 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 45 ierr = PetscLogFlops(6*m);CHKERRQ(ierr); 46 PetscFunctionReturn(0); 47 } 48 #undef __FUNCT__ 49 #define __FUNCT__ "PCApply_PBJacobi_3" 50 static PetscErrorCode PCApply_PBJacobi_3(PC pc,Vec x,Vec y) 51 { 52 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 53 PetscErrorCode ierr; 54 PetscInt i,m = jac->mbs; 55 PetscScalar *diag = jac->diag,x0,x1,x2,*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[3*i]; x1 = xx[3*i+1]; x2 = xx[3*i+2]; 62 yy[3*i] = diag[0]*x0 + diag[3]*x1 + diag[6]*x2; 63 yy[3*i+1] = diag[1]*x0 + diag[4]*x1 + diag[7]*x2; 64 yy[3*i+2] = diag[2]*x0 + diag[5]*x1 + diag[8]*x2; 65 diag += 9; 66 } 67 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 68 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 69 ierr = PetscLogFlops(15*m);CHKERRQ(ierr); 70 PetscFunctionReturn(0); 71 } 72 #undef __FUNCT__ 73 #define __FUNCT__ "PCApply_PBJacobi_4" 74 static PetscErrorCode PCApply_PBJacobi_4(PC pc,Vec x,Vec y) 75 { 76 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 77 PetscErrorCode ierr; 78 PetscInt i,m = jac->mbs; 79 PetscScalar *diag = jac->diag,x0,x1,x2,x3,*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[4*i]; x1 = xx[4*i+1]; x2 = xx[4*i+2]; x3 = xx[4*i+3]; 86 yy[4*i] = diag[0]*x0 + diag[4]*x1 + diag[8]*x2 + diag[12]*x3; 87 yy[4*i+1] = diag[1]*x0 + diag[5]*x1 + diag[9]*x2 + diag[13]*x3; 88 yy[4*i+2] = diag[2]*x0 + diag[6]*x1 + diag[10]*x2 + diag[14]*x3; 89 yy[4*i+3] = diag[3]*x0 + diag[7]*x1 + diag[11]*x2 + diag[15]*x3; 90 diag += 16; 91 } 92 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 93 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 94 ierr = PetscLogFlops(28*m);CHKERRQ(ierr); 95 PetscFunctionReturn(0); 96 } 97 #undef __FUNCT__ 98 #define __FUNCT__ "PCApply_PBJacobi_5" 99 static PetscErrorCode PCApply_PBJacobi_5(PC pc,Vec x,Vec y) 100 { 101 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 102 PetscErrorCode ierr; 103 PetscInt i,m = jac->mbs; 104 PetscScalar *diag = jac->diag,x0,x1,x2,x3,x4,*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[5*i]; x1 = xx[5*i+1]; x2 = xx[5*i+2]; x3 = xx[5*i+3]; x4 = xx[5*i+4]; 111 yy[5*i] = diag[0]*x0 + diag[5]*x1 + diag[10]*x2 + diag[15]*x3 + diag[20]*x4; 112 yy[5*i+1] = diag[1]*x0 + diag[6]*x1 + diag[11]*x2 + diag[16]*x3 + diag[21]*x4; 113 yy[5*i+2] = diag[2]*x0 + diag[7]*x1 + diag[12]*x2 + diag[17]*x3 + diag[22]*x4; 114 yy[5*i+3] = diag[3]*x0 + diag[8]*x1 + diag[13]*x2 + diag[18]*x3 + diag[23]*x4; 115 yy[5*i+4] = diag[4]*x0 + diag[9]*x1 + diag[14]*x2 + diag[19]*x3 + diag[24]*x4; 116 diag += 25; 117 } 118 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 119 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 120 ierr = PetscLogFlops(45*m);CHKERRQ(ierr); 121 PetscFunctionReturn(0); 122 } 123 /* -------------------------------------------------------------------------- */ 124 #undef __FUNCT__ 125 #define __FUNCT__ "PCSetUp_PBJacobi" 126 static PetscErrorCode PCSetUp_PBJacobi(PC pc) 127 { 128 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 129 PetscErrorCode ierr; 130 PetscMPIInt size; 131 PetscTruth seqbaij,mpibaij,baij; 132 Mat A = pc->pmat; 133 Mat_SeqBAIJ *a; 134 135 PetscFunctionBegin; 136 ierr = PetscTypeCompare((PetscObject)pc->pmat,MATSEQBAIJ,&seqbaij);CHKERRQ(ierr); 137 ierr = PetscTypeCompare((PetscObject)pc->pmat,MATMPIBAIJ,&mpibaij);CHKERRQ(ierr); 138 ierr = PetscTypeCompare((PetscObject)pc->pmat,MATBAIJ,&baij);CHKERRQ(ierr); 139 if (!seqbaij && !mpibaij && !baij) { 140 SETERRQ(PETSC_ERR_SUP,"Currently only supports BAIJ matrices"); 141 } 142 ierr = MPI_Comm_size(pc->comm,&size);CHKERRQ(ierr); 143 if (mpibaij || (baij && (size > 1))) A = ((Mat_MPIBAIJ*)A->data)->A; 144 if (A->rmap.n != A->cmap.n) SETERRQ(PETSC_ERR_SUP,"Supported only for square matrices and square storage"); 145 146 ierr = MatSeqBAIJInvertBlockDiagonal(A);CHKERRQ(ierr); 147 a = (Mat_SeqBAIJ*)A->data; 148 jac->diag = a->idiag; 149 jac->bs = A->rmap.bs; 150 jac->mbs = a->mbs; 151 switch (jac->bs){ 152 case 2: 153 pc->ops->apply = PCApply_PBJacobi_2; 154 break; 155 case 3: 156 pc->ops->apply = PCApply_PBJacobi_3; 157 break; 158 case 4: 159 pc->ops->apply = PCApply_PBJacobi_4; 160 break; 161 case 5: 162 pc->ops->apply = PCApply_PBJacobi_5; 163 break; 164 default: 165 SETERRQ1(PETSC_ERR_SUP,"not supported for block size %D",jac->bs); 166 } 167 168 PetscFunctionReturn(0); 169 } 170 /* -------------------------------------------------------------------------- */ 171 #undef __FUNCT__ 172 #define __FUNCT__ "PCDestroy_PBJacobi" 173 static PetscErrorCode PCDestroy_PBJacobi(PC pc) 174 { 175 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 176 PetscErrorCode ierr; 177 178 PetscFunctionBegin; 179 /* 180 Free the private data structure that was hanging off the PC 181 */ 182 ierr = PetscFree(jac);CHKERRQ(ierr); 183 PetscFunctionReturn(0); 184 } 185 /* -------------------------------------------------------------------------- */ 186 /*MC 187 PCPBJACOBI - Point block Jacobi 188 189 Level: beginner 190 191 Concepts: point block Jacobi 192 193 Notes: Only implemented for the BAIJ matrix formats. 194 195 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC 196 197 M*/ 198 199 EXTERN_C_BEGIN 200 #undef __FUNCT__ 201 #define __FUNCT__ "PCCreate_PBJacobi" 202 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_PBJacobi(PC pc) 203 { 204 PC_PBJacobi *jac; 205 PetscErrorCode ierr; 206 207 PetscFunctionBegin; 208 209 /* 210 Creates the private data structure for this preconditioner and 211 attach it to the PC object. 212 */ 213 ierr = PetscNewLog(pc,PC_PBJacobi,&jac);CHKERRQ(ierr); 214 pc->data = (void*)jac; 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