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_2" 26 static PetscErrorCode PCApply_PBJacobi_2(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 PetscScalar 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.0*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 const MatScalar *diag = jac->diag; 56 PetscScalar x0,x1,x2,*xx,*yy; 57 58 PetscFunctionBegin; 59 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 60 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 61 for (i=0; i<m; i++) { 62 x0 = xx[3*i]; x1 = xx[3*i+1]; x2 = xx[3*i+2]; 63 yy[3*i] = diag[0]*x0 + diag[3]*x1 + diag[6]*x2; 64 yy[3*i+1] = diag[1]*x0 + diag[4]*x1 + diag[7]*x2; 65 yy[3*i+2] = diag[2]*x0 + diag[5]*x1 + diag[8]*x2; 66 diag += 9; 67 } 68 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 69 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 70 ierr = PetscLogFlops(15.0*m);CHKERRQ(ierr); 71 PetscFunctionReturn(0); 72 } 73 #undef __FUNCT__ 74 #define __FUNCT__ "PCApply_PBJacobi_4" 75 static PetscErrorCode PCApply_PBJacobi_4(PC pc,Vec x,Vec y) 76 { 77 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 78 PetscErrorCode ierr; 79 PetscInt i,m = jac->mbs; 80 const MatScalar *diag = jac->diag; 81 PetscScalar x0,x1,x2,x3,*xx,*yy; 82 83 PetscFunctionBegin; 84 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 85 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 86 for (i=0; i<m; i++) { 87 x0 = xx[4*i]; x1 = xx[4*i+1]; x2 = xx[4*i+2]; x3 = xx[4*i+3]; 88 yy[4*i] = diag[0]*x0 + diag[4]*x1 + diag[8]*x2 + diag[12]*x3; 89 yy[4*i+1] = diag[1]*x0 + diag[5]*x1 + diag[9]*x2 + diag[13]*x3; 90 yy[4*i+2] = diag[2]*x0 + diag[6]*x1 + diag[10]*x2 + diag[14]*x3; 91 yy[4*i+3] = diag[3]*x0 + diag[7]*x1 + diag[11]*x2 + diag[15]*x3; 92 diag += 16; 93 } 94 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 95 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 96 ierr = PetscLogFlops(28.0*m);CHKERRQ(ierr); 97 PetscFunctionReturn(0); 98 } 99 #undef __FUNCT__ 100 #define __FUNCT__ "PCApply_PBJacobi_5" 101 static PetscErrorCode PCApply_PBJacobi_5(PC pc,Vec x,Vec y) 102 { 103 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 104 PetscErrorCode ierr; 105 PetscInt i,m = jac->mbs; 106 const MatScalar *diag = jac->diag; 107 PetscScalar x0,x1,x2,x3,x4,*xx,*yy; 108 109 PetscFunctionBegin; 110 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 111 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 112 for (i=0; i<m; i++) { 113 x0 = xx[5*i]; x1 = xx[5*i+1]; x2 = xx[5*i+2]; x3 = xx[5*i+3]; x4 = xx[5*i+4]; 114 yy[5*i] = diag[0]*x0 + diag[5]*x1 + diag[10]*x2 + diag[15]*x3 + diag[20]*x4; 115 yy[5*i+1] = diag[1]*x0 + diag[6]*x1 + diag[11]*x2 + diag[16]*x3 + diag[21]*x4; 116 yy[5*i+2] = diag[2]*x0 + diag[7]*x1 + diag[12]*x2 + diag[17]*x3 + diag[22]*x4; 117 yy[5*i+3] = diag[3]*x0 + diag[8]*x1 + diag[13]*x2 + diag[18]*x3 + diag[23]*x4; 118 yy[5*i+4] = diag[4]*x0 + diag[9]*x1 + diag[14]*x2 + diag[19]*x3 + diag[24]*x4; 119 diag += 25; 120 } 121 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 122 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 123 ierr = PetscLogFlops(45.0*m);CHKERRQ(ierr); 124 PetscFunctionReturn(0); 125 } 126 /* -------------------------------------------------------------------------- */ 127 #undef __FUNCT__ 128 #define __FUNCT__ "PCSetUp_PBJacobi" 129 static PetscErrorCode PCSetUp_PBJacobi(PC pc) 130 { 131 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 132 PetscErrorCode ierr; 133 PetscMPIInt size; 134 PetscBool seqbaij,mpibaij,baij; 135 Mat A = pc->pmat; 136 Mat_SeqBAIJ *a; 137 138 PetscFunctionBegin; 139 ierr = PetscTypeCompare((PetscObject)pc->pmat,MATSEQBAIJ,&seqbaij);CHKERRQ(ierr); 140 ierr = PetscTypeCompare((PetscObject)pc->pmat,MATMPIBAIJ,&mpibaij);CHKERRQ(ierr); 141 ierr = PetscTypeCompare((PetscObject)pc->pmat,MATBAIJ,&baij);CHKERRQ(ierr); 142 if (!seqbaij && !mpibaij && !baij) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Currently only supports BAIJ matrices"); 143 ierr = MPI_Comm_size(((PetscObject)pc)->comm,&size);CHKERRQ(ierr); 144 if (mpibaij || (baij && (size > 1))) A = ((Mat_MPIBAIJ*)A->data)->A; 145 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Supported only for square matrices and square storage"); 146 147 ierr = MatSeqBAIJInvertBlockDiagonal(A);CHKERRQ(ierr); 148 a = (Mat_SeqBAIJ*)A->data; 149 jac->diag = a->idiag; 150 jac->bs = A->rmap->bs; 151 jac->mbs = a->mbs; 152 switch (jac->bs){ 153 case 2: 154 pc->ops->apply = PCApply_PBJacobi_2; 155 break; 156 case 3: 157 pc->ops->apply = PCApply_PBJacobi_3; 158 break; 159 case 4: 160 pc->ops->apply = PCApply_PBJacobi_4; 161 break; 162 case 5: 163 pc->ops->apply = PCApply_PBJacobi_5; 164 break; 165 default: 166 SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"not supported for block size %D",jac->bs); 167 } 168 169 PetscFunctionReturn(0); 170 } 171 /* -------------------------------------------------------------------------- */ 172 #undef __FUNCT__ 173 #define __FUNCT__ "PCDestroy_PBJacobi" 174 static PetscErrorCode PCDestroy_PBJacobi(PC pc) 175 { 176 PetscErrorCode ierr; 177 178 PetscFunctionBegin; 179 /* 180 Free the private data structure that was hanging off the PC 181 */ 182 ierr = PetscFree(pc->data);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 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