[petsc-users] MATBDIAG in > 3.1

[Satish Balay]

On Tue, 6 Dec 2011, Matthew Knepley wrote:

> On Tue, Dec 6, 2011 at 7:43 AM, Parker, Andrew (UK Filton) <
> Andrew.Parker2 at baesystems.com> wrote:
> 
> > So having got petsc 3.2 install and all its deps, I now get this….****
> >
> > ** **
> >
> > [0]PETSC ERROR: PCSetUp_PBJacobi() line 207 in
> > src/ksp/pc/impls/pbjacobi/pbjacobi.c not supported for block size 7****
> >
> > ** **
> >
> > ** **
> >
> > Any thoughts?
> >
> 
> Sorry about that. There was confusion on our end. I have added the code you
> need to 3.2. It will
> go out with the next patch release. Alternatively, you can retrieve the
> source tarball from the
> repository site or use the repository directly.


I'm attaching this change as a patchfile - which can be applied to the current
petsc-3.2 sources you have.

cd petsc-3.2
patch -Np1 < bpjacobi.patch

Satish
-------------- next part --------------
tree baf1a1164d63
parent d7b881216ffe
author Matthew G Knepley <knepley at gmail.com> 1323197905 21600
committer Matthew G Knepley <knepley at gmail.com> 1323197905 21600
revision 20397
branch default

Added support for block size 7 to PBJACOBI
diff --git a/src/ksp/pc/impls/pbjacobi/pbjacobi.c b/src/ksp/pc/impls/pbjacobi/pbjacobi.c
--- a/src/ksp/pc/impls/pbjacobi/pbjacobi.c
+++ b/src/ksp/pc/impls/pbjacobi/pbjacobi.c
@@ -169,6 +169,35 @@
   ierr = PetscLogFlops(66.0*m);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
+#undef __FUNCT__  
+#define __FUNCT__ "PCApply_PBJacobi_7"
+static PetscErrorCode PCApply_PBJacobi_7(PC pc,Vec x,Vec y)
+{
+  PC_PBJacobi     *jac = (PC_PBJacobi*)pc->data;
+  PetscErrorCode  ierr;
+  PetscInt        i,m = jac->mbs;
+  const MatScalar *diag = jac->diag;
+  PetscScalar     x0,x1,x2,x3,x4,x5,x6,*xx,*yy;
+
+  PetscFunctionBegin;
+  ierr = VecGetArray(x,&xx);CHKERRQ(ierr);
+  ierr = VecGetArray(y,&yy);CHKERRQ(ierr);
+  for (i=0; i<m; i++) {
+    x0 = xx[7*i]; x1 = xx[7*i+1]; x2 = xx[7*i+2]; x3 = xx[7*i+3]; x4 = xx[7*i+4]; x5 = xx[7*i+5]; x6 = xx[7*i+6];
+    yy[7*i]   = diag[0]*x0 + diag[7]*x1  + diag[14]*x2  + diag[21]*x3 + diag[28]*x4 + diag[35]*x5 + diag[42]*x6;
+    yy[7*i+1] = diag[1]*x0 + diag[8]*x1  + diag[15]*x2  + diag[22]*x3 + diag[29]*x4 + diag[36]*x5 + diag[43]*x6;
+    yy[7*i+2] = diag[2]*x0 + diag[9]*x1  + diag[16]*x2  + diag[23]*x3 + diag[30]*x4 + diag[37]*x5 + diag[44]*x6;
+    yy[7*i+3] = diag[3]*x0 + diag[10]*x1 + diag[17]*x2  + diag[24]*x3 + diag[31]*x4 + diag[38]*x5 + diag[45]*x6;
+    yy[7*i+4] = diag[4]*x0 + diag[11]*x1 + diag[18]*x2  + diag[25]*x3 + diag[32]*x4 + diag[39]*x5 + diag[46]*x6;
+    yy[7*i+5] = diag[5]*x0 + diag[12]*x1 + diag[19]*x2  + diag[26]*x3 + diag[33]*x4 + diag[40]*x5 + diag[47]*x6;
+    yy[7*i+6] = diag[6]*x0 + diag[13]*x1 + diag[20]*x2  + diag[27]*x3 + diag[34]*x4 + diag[41]*x5 + diag[48]*x6;
+    diag     += 49;
+  }
+  ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr);
+  ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr);
+  ierr = PetscLogFlops(80.0*m);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
 /* -------------------------------------------------------------------------- */
 #undef __FUNCT__  
 #define __FUNCT__ "PCSetUp_PBJacobi"
@@ -203,6 +232,9 @@
     case 6:
       pc->ops->apply = PCApply_PBJacobi_6;
       break;
+    case 7:
+      pc->ops->apply = PCApply_PBJacobi_7;
+      break;
     default: 
       SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"not supported for block size %D",jac->bs);
   }