1. Vijay Mahadevan
  2. petsc

Commits

BarryFSmith  committed 188559c Merge

Merge branch 'barry/add-vecstrideblock'

  • Participants
  • Parent commits 0157cb0, bdeb6c8
  • Branches master

Comments (0)

Files changed (8)

File include/petsc-private/vecimpl.h

View file
  • Ignore whitespace
   PetscErrorCode (*restoresubvector)(Vec,IS,Vec*);
   PetscErrorCode (*getarrayread)(Vec,const PetscScalar**);
   PetscErrorCode (*restorearrayread)(Vec,const PetscScalar**);
+  PetscErrorCode (*stridesubsetgather)(Vec,PetscInt,const PetscInt[],const PetscInt[],Vec,InsertMode);
+  PetscErrorCode (*stridesubsetscatter)(Vec,PetscInt,const PetscInt[],const PetscInt[],Vec,InsertMode);
 };
 
 /*
 PETSC_INTERN PetscErrorCode VecStrideGather_Default(Vec,PetscInt,Vec,InsertMode);
 PETSC_INTERN PetscErrorCode VecStrideScatter_Default(Vec,PetscInt,Vec,InsertMode);
 PETSC_INTERN PetscErrorCode VecReciprocal_Default(Vec);
+PETSC_INTERN PetscErrorCode VecStrideSubSetGather_Default(Vec,PetscInt,const PetscInt[],const PetscInt[],Vec,InsertMode);
+PETSC_INTERN PetscErrorCode VecStrideSubSetScatter_Default(Vec,PetscInt,const PetscInt[],const PetscInt[],Vec,InsertMode);
 
 #if defined(PETSC_HAVE_MATLAB_ENGINE)
 PETSC_EXTERN PetscErrorCode VecMatlabEnginePut_Default(PetscObject,void*);

File include/petscvec.h

View file
  • Ignore whitespace
 PETSC_EXTERN PetscErrorCode VecStrideGatherAll(Vec,Vec[],InsertMode);
 PETSC_EXTERN PetscErrorCode VecStrideScatterAll(Vec[],Vec,InsertMode);
 
+PETSC_EXTERN PetscErrorCode VecStrideSubSetScatter(Vec,PetscInt,const PetscInt[],const PetscInt[],Vec,InsertMode);
+PETSC_EXTERN PetscErrorCode VecStrideSubSetGather(Vec,PetscInt,const PetscInt[],const PetscInt[],Vec,InsertMode);
+
 PETSC_EXTERN PetscErrorCode VecSetValues(Vec,PetscInt,const PetscInt[],const PetscScalar[],InsertMode);
 PETSC_EXTERN PetscErrorCode VecGetValues(Vec,PetscInt,const PetscInt[],PetscScalar[]);
 PETSC_EXTERN PetscErrorCode VecAssemblyBegin(Vec);

File src/vec/vec/examples/tests/ex45.c

View file
  • Ignore whitespace
+
+static char help[] = "Demonstrates VecStrideSubSetScatter() and VecStrideSubSetGather().\n\n";
+
+/*T
+   Concepts: vectors^sub-vectors;
+   Processors: n
+
+   Allows one to easily pull out some components of a multi-component vector and put them in another vector.
+
+   Note that these are special cases of VecScatter
+T*/
+
+/*
+  Include "petscvec.h" so that we can use vectors.  Note that this file
+  automatically includes:
+     petscsys.h       - base PETSc routines   petscis.h     - index sets
+     petscviewer.h - viewers
+*/
+
+#include <petscvec.h>
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc,char **argv)
+{
+  Vec            v,s;
+  PetscInt       i,start,end,n = 8;
+  PetscErrorCode ierr;
+  PetscScalar    value;
+  const PetscInt vidx[] = {1,2},sidx[] = {1,0};
+
+  PetscInitialize(&argc,&argv,(char*)0,help);
+  ierr = PetscOptionsGetInt(NULL,"-n",&n,NULL);CHKERRQ(ierr);
+
+  /*
+      Create multi-component vector with 4 components
+  */
+  ierr = VecCreate(PETSC_COMM_WORLD,&v);CHKERRQ(ierr);
+  ierr = VecSetSizes(v,PETSC_DECIDE,n);CHKERRQ(ierr);
+  ierr = VecSetBlockSize(v,4);CHKERRQ(ierr);
+  ierr = VecSetFromOptions(v);CHKERRQ(ierr);
+
+  /*
+      Create double-component vectors
+  */
+  ierr = VecCreate(PETSC_COMM_WORLD,&s);CHKERRQ(ierr);
+  ierr = VecSetSizes(s,PETSC_DECIDE,n/2);CHKERRQ(ierr);
+  ierr = VecSetBlockSize(s,2);CHKERRQ(ierr);
+  ierr = VecSetFromOptions(s);CHKERRQ(ierr);
+
+  /*
+     Set the vector values
+  */
+  ierr = VecGetOwnershipRange(v,&start,&end);CHKERRQ(ierr);
+  for (i=start; i<end; i++) {
+    value = i;
+    ierr  = VecSetValues(v,1,&i,&value,INSERT_VALUES);CHKERRQ(ierr);
+  }
+
+  /*
+     Get the components from the large multi-component vector to the small multi-component vector,
+     scale the smaller vector and then move values back to the large vector
+  */
+  ierr = VecStrideSubSetGather(v,PETSC_DETERMINE,vidx,NULL,s,INSERT_VALUES);CHKERRQ(ierr);
+  ierr = VecView(s,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
+  ierr = VecScale(s,100.0);CHKERRQ(ierr);
+
+  ierr = VecStrideSubSetScatter(s,PETSC_DETERMINE,NULL,vidx,v,ADD_VALUES);CHKERRQ(ierr);
+  ierr = VecView(v,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
+
+  /*
+     Get the components from the large multi-component vector to the small multi-component vector,
+     scale the smaller vector and then move values back to the large vector
+  */
+  ierr = VecStrideSubSetGather(v,2,vidx,sidx,s,INSERT_VALUES);CHKERRQ(ierr);
+  ierr = VecView(s,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
+  ierr = VecScale(s,100.0);CHKERRQ(ierr);
+
+  ierr = VecStrideSubSetScatter(s,2,sidx,vidx,v,ADD_VALUES);CHKERRQ(ierr);
+  ierr = VecView(v,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
+
+  /*
+     Free work space.  All PETSc objects should be destroyed when they
+     are no longer needed.
+  */
+  ierr = VecDestroy(&v);CHKERRQ(ierr);
+  ierr = VecDestroy(&s);CHKERRQ(ierr);
+  ierr = PetscFinalize();
+  return 0;
+}
+

File src/vec/vec/examples/tests/makefile

View file
  • Ignore whitespace
 EXAMPLESC       = ex1.c ex2.c ex3.c ex4.c ex5.c ex6.c ex7.c ex8.c ex9.c ex10.c \
                 ex11.c ex12.c ex14.c ex15.c ex16.c ex17.c ex18.c ex21.c ex22.c \
                 ex23.c ex24.c ex25.c ex28.c ex29.c ex31.c ex33.c ex34.c ex35.c \
-                ex36.c ex37.c ex38.c ex39.c ex40.c ex41.c ex42.c
+                ex36.c ex37.c ex38.c ex39.c ex40.c ex41.c ex42.c ex45.c
 EXAMPLESF       = ex17f.F ex19f.F ex20f.F ex30f.F ex32f.F
 MANSEC          = Vec
 
 	-${CLINKER} -o ex44 ex44.o ${PETSC_VEC_LIB}
 	${RM} -f ex44.o
 
+ex45: ex45.o  chkopts
+	-${CLINKER} -o ex45 ex45.o ${PETSC_VEC_LIB}
+	${RM} -f ex45.o
+
 #--------------------------------------------------------------------------
 runex1:
 	-@${MPIEXEC} -n 1 ./ex1 > ex1_1.tmp 2>&1;\
 	   ${DIFF} output/ex44.out ex44.tmp || printf "${PWD}\nPossible problem with ex44, diffs above\n=========================================\n"; \
 	   ${RM} -f ex44.tmp
 
+runex45:
+	-@${MPIEXEC} -n 2 ./ex45  > ex45.tmp 2>&1;\
+	   ${DIFF} output/ex45_1.out ex45.tmp || printf "${PWD}\nPossible problem with ex45, diffs above\n=========================================\n"; \
+	   ${RM} -f ex45.tmp
+
 TESTEXAMPLES_C		    = ex1.PETSc runex1 ex1.rm ex2.PETSc runex2 ex2.rm ex3.PETSc runex3 runex3_2 ex3.rm \
                               ex4.PETSc runex4 ex4.rm ex5.PETSc ex5.rm ex6.PETSc runex6 ex6.rm ex7.PETSc \
                               runex7 ex7.rm ex8.PETSc runex8 ex8.rm ex9.PETSc runex9 ex9.rm ex11.PETSc runex11 \
                               ex14.rm ex15.PETSc runex15 ex15.rm ex16.PETSc runex16 ex16.rm ex17.PETSc runex17 \
                               ex17.rm ex21.PETSc runex21 runex21_2 ex21.rm ex25.PETSc runex25 ex25.rm ex29.PETSc \
                               runex29 ex29.rm ex34.PETSc runex34 ex34.rm ex36.PETSc runex36 ex36.rm \
-                              ex37.PETSc runex37 runex37_2 runex37_3 runex37_4  ex37.rm ex38.PETSc runex38 ex38.rm
+                              ex37.PETSc runex37 runex37_2 runex37_3 runex37_4  ex37.rm ex38.PETSc runex38 ex38.rm ex45.PETSc runex45 ex45.rm
 TESTEXAMPLES_C_X	    = ex10.PETSc runex10 ex10.rm ex22.PETSc runex22 ex22.rm ex23.PETSc runex23 ex23.rm \
                               ex24.PETSc runex24 ex24.rm ex28.PETSc runex28 runex28_2 ex28.rm ex33.PETSc runex33 ex33.rm
 TESTEXAMPLES_FORTRAN	    = ex17f.PETSc runex17f ex17f.rm ex19f.PETSc ex19f.rm ex20f.PETSc ex20f.rm ex30f.PETSc \

File src/vec/vec/examples/tests/output/ex45_1.out

View file
  • Ignore whitespace
+Vec Object: 2 MPI processes
+  type: mpi
+Process [0]
+1
+2
+Process [1]
+5
+6
+Vec Object: 2 MPI processes
+  type: mpi
+Process [0]
+0
+101
+202
+3
+Process [1]
+4
+505
+606
+7
+Vec Object: 2 MPI processes
+  type: mpi
+Process [0]
+202
+101
+Process [1]
+606
+505
+Vec Object: 2 MPI processes
+  type: mpi
+Process [0]
+0
+10201
+20402
+3
+Process [1]
+4
+51005
+61206
+7

File src/vec/vec/impls/mpi/pbvec.c

View file
  • Ignore whitespace
                                 0,
                                 0,
                                 0,
-                                0
+                                0,
+                                VecStrideSubSetGather_Default,
+                                VecStrideSubSetScatter_Default
 };
 
 #undef __FUNCT__

File src/vec/vec/impls/seq/bvec2.c

View file
  • Ignore whitespace
                                0,
                                0,
                                0,
-                               0
+                               0,
+                               VecStrideSubSetGather_Default,
+                               VecStrideSubSetScatter_Default
 };
 
 

File src/vec/vec/utils/vinv.c

View file
  • Ignore whitespace
    Concepts: scatter^into strided vector
 
 .seealso: VecStrideNorm(), VecStrideGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
-          VecStrideScatterAll()
+          VecStrideScatterAll(), VecStrideSubSetScatter(), VecStrideSubSetGather()
 @*/
 PetscErrorCode  VecStrideScatter(Vec s,PetscInt start,Vec v,InsertMode addv)
 {
 }
 
 #undef __FUNCT__
+#define __FUNCT__ "VecStrideSubSetGather"
+/*@
+   VecStrideSubSetGather - Gathers a subset of components from a multi-component vector into
+   another vector.
+
+   Collective on Vec
+
+   Input Parameter:
++  v - the vector
+.  nidx - the number of indices
+.  idxv - the indices of the components 0 <= idxv[0] ...idxv[nidx-1] < bs(v), they need not be sorted
+.  idxs - the indices of the components 0 <= idxs[0] ...idxs[nidx-1] < bs(s), they need not be sorted, may be null if nidx == bs(s) or is PETSC_DETERMINE
+-  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
+
+   Output Parameter:
+.  s - the location where the subvector is stored
+
+   Notes:
+   One must call VecSetBlockSize() on both vectors before this routine to set the stride
+   information, or use a vector created from a multicomponent DMDA.
+
+
+   The parallel layout of the vector and the subvector must be the same;
+
+   Not optimized; could be easily
+
+   Level: advanced
+
+   Concepts: gather^into strided vector
+
+.seealso: VecStrideNorm(), VecStrideScatter(), VecStrideGather(), VecStrideSubSetScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
+          VecStrideScatterAll()
+@*/
+PetscErrorCode  VecStrideSubSetGather(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
+{
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(v,VEC_CLASSID,1);
+  PetscValidHeaderSpecific(s,VEC_CLASSID,5);
+  if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
+  if (!v->ops->stridesubsetgather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
+  ierr = (*v->ops->stridesubsetgather)(v,nidx,idxv,idxs,s,addv);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "VecStrideSubSetScatter"
+/*@
+   VecStrideSubSetScatter - Scatters components from a vector into a subset of components of a multi-component vector.
+
+   Collective on Vec
+
+   Input Parameter:
++  s - the smaller-component vector
+.  nidx - the number of indices in idx
+.  idxs - the indices of the components in the smaller-component vector, 0 <= idxs[0] ...idxs[nidx-1] < bs(s) they need not be sorted, may be null if nidx == bs(s) or is PETSC_DETERMINE
+.  idxv - the indices of the components in the larger-component vector, 0 <= idx[0] ...idx[nidx-1] < bs(v) they need not be sorted
+-  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
+
+   Output Parameter:
+.  v - the location where the subvector is into scattered (the multi-component vector)
+
+   Notes:
+   One must call VecSetBlockSize() on the vectors before this
+   routine to set the stride  information, or use a vector created from a multicomponent DMDA.
+
+   The parallel layout of the vector and the subvector must be the same;
+
+   Not optimized; could be easily
+
+   Level: advanced
+
+   Concepts: scatter^into strided vector
+
+.seealso: VecStrideNorm(), VecStrideGather(), VecStrideGather(), VecStrideSubSetGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
+          VecStrideScatterAll()
+@*/
+PetscErrorCode  VecStrideSubSetScatter(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
+{
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(s,VEC_CLASSID,1);
+  PetscValidHeaderSpecific(v,VEC_CLASSID,5);
+  if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
+  if (!v->ops->stridesubsetscatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
+  ierr = (*v->ops->stridesubsetscatter)(s,nidx,idxs,idxv,v,addv);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
 #define __FUNCT__ "VecStrideGather_Default"
 PetscErrorCode  VecStrideGather_Default(Vec v,PetscInt start,Vec s,InsertMode addv)
 {
   PetscErrorCode ierr;
   PetscInt       i,n,bs,ns;
-  PetscScalar    *x,*y;
+  const PetscScalar *x;
+  PetscScalar       *y;
 
   PetscFunctionBegin;
   ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
   ierr = VecGetLocalSize(s,&ns);CHKERRQ(ierr);
-  ierr = VecGetArray(v,&x);CHKERRQ(ierr);
+  ierr = VecGetArrayRead(v,&x);CHKERRQ(ierr);
   ierr = VecGetArray(s,&y);CHKERRQ(ierr);
 
   bs = v->map->bs;
 #endif
   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
 
-  ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
+  ierr = VecRestoreArrayRead(v,&x);CHKERRQ(ierr);
   ierr = VecRestoreArray(s,&y);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 #define __FUNCT__ "VecStrideScatter_Default"
 PetscErrorCode  VecStrideScatter_Default(Vec s,PetscInt start,Vec v,InsertMode addv)
 {
-  PetscErrorCode ierr;
-  PetscInt       i,n,bs,ns;
-  PetscScalar    *x,*y;
+  PetscErrorCode    ierr;
+  PetscInt          i,n,bs,ns;
+  PetscScalar       *x;
+  const PetscScalar *y;
 
   PetscFunctionBegin;
   ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
   ierr = VecGetLocalSize(s,&ns);CHKERRQ(ierr);
   ierr = VecGetArray(v,&x);CHKERRQ(ierr);
-  ierr = VecGetArray(s,&y);CHKERRQ(ierr);
+  ierr = VecGetArrayRead(s,&y);CHKERRQ(ierr);
 
   bs = v->map->bs;
   if (n != ns*bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Subvector length * blocksize %D not correct for scatter to multicomponent vector %D",ns*bs,n);
   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
 
   ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
+  ierr = VecRestoreArrayRead(s,&y);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "VecStrideSubSetGather_Default"
+PetscErrorCode  VecStrideSubSetGather_Default(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
+{
+  PetscErrorCode    ierr;
+  PetscInt          i,j,n,bs,bss,ns;
+  const PetscScalar *x;
+  PetscScalar       *y;
+
+  PetscFunctionBegin;
+  ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
+  ierr = VecGetLocalSize(s,&ns);CHKERRQ(ierr);
+  ierr = VecGetArrayRead(v,&x);CHKERRQ(ierr);
+  ierr = VecGetArray(s,&y);CHKERRQ(ierr);
+
+  bs  = v->map->bs;
+  bss = s->map->bs;
+  n  =  n/bs;
+
+#if defined(PETSC_DEBUG)
+  if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
+  for (j=0; j<nidx; j++) {
+    if (idxv[j] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idxv[j]);
+    if (idxv[j] >= bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is greater than or equal to vector blocksize %D",j,idxv[j],bs);
+  }
+  if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not gathering into all locations");
+#endif
+
+  if (addv == INSERT_VALUES) {
+    if (!idxs) {
+      for (i=0; i<n; i++) {
+        for (j=0; j<bss; j++) y[bss*i+j] = x[bs*i+idxv[j]];
+      }
+    } else {
+      for (i=0; i<n; i++) {
+        for (j=0; j<bss; j++) y[bss*i+idxs[j]] = x[bs*i+idxv[j]];
+      }
+    }
+  } else if (addv == ADD_VALUES) {
+    if (!idxs) {
+      for (i=0; i<n; i++) {
+        for (j=0; j<bss; j++) y[bss*i+j] += x[bs*i+idxv[j]];
+      }
+    } else {
+      for (i=0; i<n; i++) {
+        for (j=0; j<bss; j++) y[bss*i+idxs[j]] += x[bs*i+idxv[j]];
+      }
+    }
+#if !defined(PETSC_USE_COMPLEX)
+  } else if (addv == MAX_VALUES) {
+    if (!idxs) {
+      for (i=0; i<n; i++) {
+        for (j=0; j<bss; j++) y[bss*i+j] = PetscMax(y[bss*i+j],x[bs*i+idxv[j]]);
+      }
+    } else {
+      for (i=0; i<n; i++) {
+        for (j=0; j<bss; j++) y[bss*i+idxs[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i+idxv[j]]);
+      }
+    }
+#endif
+  } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
+
+  ierr = VecRestoreArrayRead(v,&x);CHKERRQ(ierr);
   ierr = VecRestoreArray(s,&y);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
 #undef __FUNCT__
+#define __FUNCT__ "VecStrideSubSetScatter_Default"
+PetscErrorCode  VecStrideSubSetScatter_Default(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
+{
+  PetscErrorCode    ierr;
+  PetscInt          j,i,n,bs,ns,bss;
+  PetscScalar       *x;
+  const PetscScalar *y;
+
+  PetscFunctionBegin;
+  ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
+  ierr = VecGetLocalSize(s,&ns);CHKERRQ(ierr);
+  ierr = VecGetArray(v,&x);CHKERRQ(ierr);
+  ierr = VecGetArrayRead(s,&y);CHKERRQ(ierr);
+
+  bs  = v->map->bs;
+  bss = s->map->bs;
+  n  =  n/bs;
+
+#if defined(PETSC_DEBUG)
+  if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
+  for (j=0; j<bss; j++) {
+    if (idx[j] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idx[j]);
+    if (idx[j] >= bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is greater than or equal to vector blocksize %D",j,idx[j],bs);
+  }
+  if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not scattering from all locations");
+#endif
+
+  if (addv == INSERT_VALUES) {
+    if (!idxs) {
+      for (i=0; i<n; i++) {
+        for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+j];
+      }
+    } else {
+      for (i=0; i<n; i++) {
+        for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+idxs[j]];
+      }
+    }
+  } else if (addv == ADD_VALUES) {
+    if (!idxs) {
+      for (i=0; i<n; i++) {
+        for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+j];
+      }
+    } else {
+      for (i=0; i<n; i++) {
+        for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+idxs[j]];
+      }
+    }
+#if !defined(PETSC_USE_COMPLEX)
+  } else if (addv == MAX_VALUES) {
+    if (!idxs) {
+      for (i=0; i<n; i++) {
+        for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+j],x[bs*i + idxv[j]]);
+      }
+    } else {
+      for (i=0; i<n; i++) {
+        for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i + idxv[j]]);
+      }
+    }
+#endif
+  } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
+
+  ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
+  ierr = VecRestoreArrayRead(s,&y);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
 #define __FUNCT__ "VecReciprocal_Default"
 PetscErrorCode VecReciprocal_Default(Vec v)
 {