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/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 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.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 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.0*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.0*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.0*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) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Currently only supports BAIJ matrices"); 140 ierr = MPI_Comm_size(((PetscObject)pc)->comm,&size);CHKERRQ(ierr); 141 if (mpibaij || (baij && (size > 1))) A = ((Mat_MPIBAIJ*)A->data)->A; 142 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,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->rmap->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(((PetscObject)pc)->comm,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 PETSCKSP_DLLEXPORT 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 = PetscNewLog(pc,PC_PBJacobi,&jac);CHKERRQ(ierr); 212 pc->data = (void*)jac; 213 214 /* 215 Initialize the pointers to vectors to ZERO; these will be used to store 216 diagonal entries of the matrix for fast preconditioner application. 217 */ 218 jac->diag = 0; 219 220 /* 221 Set the pointers for the functions that are provided above. 222 Now when the user-level routines (such as PCApply(), PCDestroy(), etc.) 223 are called, they will automatically call these functions. Note we 224 choose not to provide a couple of these functions since they are 225 not needed. 226 */ 227 pc->ops->apply = 0; /*set depending on the block size */ 228 pc->ops->applytranspose = 0; 229 pc->ops->setup = PCSetUp_PBJacobi; 230 pc->ops->destroy = PCDestroy_PBJacobi; 231 pc->ops->setfromoptions = 0; 232 pc->ops->view = 0; 233 pc->ops->applyrichardson = 0; 234 pc->ops->applysymmetricleft = 0; 235 pc->ops->applysymmetricright = 0; 236 PetscFunctionReturn(0); 237 } 238 EXTERN_C_END 239 240 241