1 /*$Id: pbjacobi.c,v 1.4 2001/08/07 03:03:42 balay Exp $*/ 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 "src/ksp/pc/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 int 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 int PCApply_PBJacobi_2(PC pc,Vec x,Vec y) 28 { 29 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 30 int ierr,i,m = jac->mbs; 31 PetscScalar *diag = jac->diag,x0,x1,*xx,*yy; 32 33 PetscFunctionBegin; 34 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 35 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 36 for (i=0; i<m; i++) { 37 x0 = xx[2*i]; x1 = xx[2*i+1]; 38 yy[2*i] = diag[0]*x0 + diag[2]*x1; 39 yy[2*i+1] = diag[1]*x0 + diag[3]*x1; 40 diag += 4; 41 } 42 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 43 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 44 PetscLogFlops(6*m); 45 PetscFunctionReturn(0); 46 } 47 #undef __FUNCT__ 48 #define __FUNCT__ "PCApply_PBJacobi_3" 49 static int PCApply_PBJacobi_3(PC pc,Vec x,Vec y) 50 { 51 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 52 int ierr,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 PetscLogFlops(15*m); 68 PetscFunctionReturn(0); 69 } 70 #undef __FUNCT__ 71 #define __FUNCT__ "PCApply_PBJacobi_4" 72 static int PCApply_PBJacobi_4(PC pc,Vec x,Vec y) 73 { 74 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 75 int ierr,i,m = jac->mbs; 76 PetscScalar *diag = jac->diag,x0,x1,x2,x3,*xx,*yy; 77 78 PetscFunctionBegin; 79 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 80 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 81 for (i=0; i<m; i++) { 82 x0 = xx[4*i]; x1 = xx[4*i+1]; x2 = xx[4*i+2]; x3 = xx[4*i+3]; 83 yy[4*i] = diag[0]*x0 + diag[4]*x1 + diag[8]*x2 + diag[12]*x3; 84 yy[4*i+1] = diag[1]*x0 + diag[5]*x1 + diag[9]*x2 + diag[13]*x3; 85 yy[4*i+2] = diag[2]*x0 + diag[6]*x1 + diag[10]*x2 + diag[14]*x3; 86 yy[4*i+3] = diag[3]*x0 + diag[7]*x1 + diag[11]*x2 + diag[15]*x3; 87 diag += 16; 88 } 89 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 90 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 91 PetscLogFlops(28*m); 92 PetscFunctionReturn(0); 93 } 94 #undef __FUNCT__ 95 #define __FUNCT__ "PCApply_PBJacobi_5" 96 static int PCApply_PBJacobi_5(PC pc,Vec x,Vec y) 97 { 98 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 99 int ierr,i,m = jac->mbs; 100 PetscScalar *diag = jac->diag,x0,x1,x2,x3,x4,*xx,*yy; 101 102 PetscFunctionBegin; 103 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 104 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 105 for (i=0; i<m; i++) { 106 x0 = xx[5*i]; x1 = xx[5*i+1]; x2 = xx[5*i+2]; x3 = xx[5*i+3]; x4 = xx[5*i+4]; 107 yy[5*i] = diag[0]*x0 + diag[5]*x1 + diag[10]*x2 + diag[15]*x3 + diag[20]*x4; 108 yy[5*i+1] = diag[1]*x0 + diag[6]*x1 + diag[11]*x2 + diag[16]*x3 + diag[21]*x4; 109 yy[5*i+2] = diag[2]*x0 + diag[7]*x1 + diag[12]*x2 + diag[17]*x3 + diag[22]*x4; 110 yy[5*i+3] = diag[3]*x0 + diag[8]*x1 + diag[13]*x2 + diag[18]*x3 + diag[23]*x4; 111 yy[5*i+4] = diag[4]*x0 + diag[9]*x1 + diag[14]*x2 + diag[19]*x3 + diag[24]*x4; 112 diag += 25; 113 } 114 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 115 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 116 PetscLogFlops(45*m); 117 PetscFunctionReturn(0); 118 } 119 /* -------------------------------------------------------------------------- */ 120 extern int MatInvertBlockDiagonal_SeqBAIJ(Mat); 121 #undef __FUNCT__ 122 #define __FUNCT__ "PCSetUp_PBJacobi" 123 static int PCSetUp_PBJacobi(PC pc) 124 { 125 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 126 int ierr,size; 127 PetscTruth seqbaij,mpibaij,baij; 128 Mat A = pc->pmat; 129 Mat_SeqBAIJ *a; 130 131 PetscFunctionBegin; 132 ierr = PetscTypeCompare((PetscObject)pc->pmat,MATSEQBAIJ,&seqbaij);CHKERRQ(ierr); 133 ierr = PetscTypeCompare((PetscObject)pc->pmat,MATMPIBAIJ,&mpibaij);CHKERRQ(ierr); 134 ierr = PetscTypeCompare((PetscObject)pc->pmat,MATBAIJ,&baij);CHKERRQ(ierr); 135 if (!seqbaij && !mpibaij && !baij) { 136 SETERRQ(1,"Currently only supports BAIJ matrices"); 137 } 138 ierr = MPI_Comm_size(pc->comm,&size);CHKERRQ(ierr); 139 if (mpibaij || (baij && (size > 1))) A = ((Mat_MPIBAIJ*)A->data)->A; 140 if (A->m != A->n) SETERRQ(1,"Supported only for square matrices and square storage"); 141 142 ierr = MatInvertBlockDiagonal_SeqBAIJ(A);CHKERRQ(ierr); 143 a = (Mat_SeqBAIJ*)A->data; 144 jac->diag = a->idiag; 145 jac->bs = a->bs; 146 jac->mbs = a->mbs; 147 switch (a->bs){ 148 case 2: 149 pc->ops->apply = PCApply_PBJacobi_2; 150 break; 151 case 3: 152 pc->ops->apply = PCApply_PBJacobi_3; 153 break; 154 case 4: 155 pc->ops->apply = PCApply_PBJacobi_4; 156 break; 157 case 5: 158 pc->ops->apply = PCApply_PBJacobi_5; 159 break; 160 default: 161 SETERRQ1(1,"not supported for block size %d",a->bs); 162 } 163 164 PetscFunctionReturn(0); 165 } 166 /* -------------------------------------------------------------------------- */ 167 #undef __FUNCT__ 168 #define __FUNCT__ "PCDestroy_PBJacobi" 169 static int PCDestroy_PBJacobi(PC pc) 170 { 171 PC_PBJacobi *jac = (PC_PBJacobi*)pc->data; 172 int ierr; 173 174 PetscFunctionBegin; 175 /* 176 Free the private data structure that was hanging off the PC 177 */ 178 ierr = PetscFree(jac);CHKERRQ(ierr); 179 PetscFunctionReturn(0); 180 } 181 /* -------------------------------------------------------------------------- */ 182 /*MC 183 PCPBJACOBI - Point block Jacobi 184 185 Level: beginner 186 187 Concepts: point block Jacobi 188 189 Notes: Only implemented for the BAIJ matrix formats. 190 191 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC 192 193 M*/ 194 195 EXTERN_C_BEGIN 196 #undef __FUNCT__ 197 #define __FUNCT__ "PCCreate_PBJacobi" 198 int PCCreate_PBJacobi(PC pc) 199 { 200 PC_PBJacobi *jac; 201 int ierr; 202 203 PetscFunctionBegin; 204 205 /* 206 Creates the private data structure for this preconditioner and 207 attach it to the PC object. 208 */ 209 ierr = PetscNew(PC_PBJacobi,&jac);CHKERRQ(ierr); 210 pc->data = (void*)jac; 211 212 /* 213 Logs the memory usage; this is not needed but allows PETSc to 214 monitor how much memory is being used for various purposes. 215 */ 216 PetscLogObjectMemory(pc,sizeof(PC_PBJacobi)); 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