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