static PetscErrorCode PetscDrawGetSingleton_Image(PetscDraw draw,PetscDraw *sdraw)
{
  PetscImage     pimg = (PetscImage)draw->data;
  PetscImage     simg;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = PetscDrawCreate(PETSC_COMM_SELF,NULL,NULL,0,0,draw->w,draw->h,sdraw);CHKERRQ(ierr);
  ierr = PetscDrawSetType(*sdraw,PETSC_DRAW_IMAGE);CHKERRQ(ierr);
  (*sdraw)->ops->resizewindow = NULL;
  simg = (PetscImage)(*sdraw)->data;
  ierr = PetscMemcpy(simg->buffer,pimg->buffer,(size_t)(pimg->w*pimg->h));CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*
   PetscSplitReductionExtend - Double the amount of space (slots) allocated for a split reduction object.
*/
PetscErrorCode  PetscSplitReductionExtend(PetscSplitReduction *sr)
{
  PetscErrorCode ierr;
  PetscInt       maxops = sr->maxops,*reducetype = sr->reducetype;
  PetscScalar    *lvalues = sr->lvalues,*gvalues = sr->gvalues;
  void           *invecs = sr->invecs;

  PetscFunctionBegin;
  sr->maxops     = 2*maxops;
  ierr = PetscMalloc(2*2*maxops*sizeof(PetscScalar),&sr->lvalues);CHKERRQ(ierr);
  ierr = PetscMalloc(2*2*maxops*sizeof(PetscScalar),&sr->gvalues);CHKERRQ(ierr);
  ierr = PetscMalloc(2*maxops*sizeof(PetscInt),&sr->reducetype);CHKERRQ(ierr);
  ierr = PetscMalloc(2*maxops*sizeof(void*),&sr->invecs);CHKERRQ(ierr);
  ierr = PetscMemcpy(sr->lvalues,lvalues,maxops*sizeof(PetscScalar));CHKERRQ(ierr);
  ierr = PetscMemcpy(sr->gvalues,gvalues,maxops*sizeof(PetscScalar));CHKERRQ(ierr);
  ierr = PetscMemcpy(sr->reducetype,reducetype,maxops*sizeof(PetscInt));CHKERRQ(ierr);
  ierr = PetscMemcpy(sr->invecs,invecs,maxops*sizeof(void*));CHKERRQ(ierr);
  ierr = PetscFree(lvalues);CHKERRQ(ierr);
  ierr = PetscFree(gvalues);CHKERRQ(ierr);
  ierr = PetscFree(reducetype);CHKERRQ(ierr);
  ierr = PetscFree(invecs);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*@
    PetscViewerDrawSetBounds - sets the upper and lower bounds to be used in plotting

    Collective on PetscViewer

    Input Parameters:
+   viewer - the PetscViewer (created with PetscViewerDrawOpen())
.   nbounds - number of plots that can be made with this viewer, for example the dof passed to DMDACreate()
-   bounds - the actual bounds, the size of this is 2*nbounds, the values are stored in the order min F_0, max F_0, min F_1, max F_1, .....


    Options Database:
.   -draw_bounds  minF0,maxF0,minF1,maxF1

    Level: intermediate

    Notes: this determines the colors used in 2d contour plots generated with VecView() for DMDA in 2d. Any values in the vector below or above the
      bounds are moved to the bound value before plotting. In this way the color index from color to physical value remains the same for all plots generated with
      this viewer. Otherwise the color to physical value meaning changes with each new image if this is not set.

   Concepts: drawing^accessing PetscDraw context from PetscViewer
   Concepts: graphics

.seealso: PetscViewerDrawGetLG(), PetscViewerDrawGetAxis(), PetscViewerDrawOpen()
@*/
PetscErrorCode  PetscViewerDrawSetBounds(PetscViewer viewer,PetscInt nbounds,const PetscReal *bounds)
{
  PetscViewer_Draw *vdraw = (PetscViewer_Draw*)viewer->data;
  PetscErrorCode   ierr;


  PetscFunctionBegin;
  PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,1);
  vdraw->nbounds = nbounds;

  ierr = PetscMalloc(2*nbounds*sizeof(PetscReal),&vdraw->bounds);CHKERRQ(ierr);
  ierr = PetscMemcpy(vdraw->bounds,bounds,2*nbounds*sizeof(PetscReal));CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

static PetscErrorCode  ISBlockSetIndices_Block(IS is,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
{
  PetscErrorCode ierr;
  PetscInt       i,min,max;
  IS_Block       *sub = (IS_Block*)is->data;

  PetscFunctionBegin;
  if (bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"block size < 1");
  if (n < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"length < 0");
  if (n) PetscValidIntPointer(idx,3);

  ierr = PetscLayoutSetLocalSize(is->map, n*bs);CHKERRQ(ierr);
  ierr = PetscLayoutSetBlockSize(is->map, bs);CHKERRQ(ierr);
  ierr = PetscLayoutSetUp(is->map);CHKERRQ(ierr);

  if (sub->allocated) {ierr = PetscFree(sub->idx);CHKERRQ(ierr);}
  if (mode == PETSC_COPY_VALUES) {
    ierr = PetscMalloc1(n,&sub->idx);CHKERRQ(ierr);
    ierr = PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));CHKERRQ(ierr);
    ierr = PetscMemcpy(sub->idx,idx,n*sizeof(PetscInt));CHKERRQ(ierr);
    sub->allocated = PETSC_TRUE;
  } else if (mode == PETSC_OWN_POINTER) {
    sub->idx = (PetscInt*) idx;
    ierr = PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));CHKERRQ(ierr);
    sub->allocated = PETSC_TRUE;
  } else if (mode == PETSC_USE_POINTER) {
    sub->idx = (PetscInt*) idx;
    sub->allocated = PETSC_FALSE;
  }

  sub->sorted = PETSC_TRUE;
  for (i=1; i<n; i++) {
    if (idx[i] < idx[i-1]) {sub->sorted = PETSC_FALSE; break;}
  }
  if (n) {
    min = max = idx[0];
    for (i=1; i<n; i++) {
      if (idx[i] < min) min = idx[i];
      if (idx[i] > max) max = idx[i];
    }
    is->min = bs*min;
    is->max = bs*max+bs-1;
  } else {
    is->min = PETSC_MAX_INT;
    is->max = PETSC_MIN_INT;
  }
  is->isperm     = PETSC_FALSE;
  is->isidentity = PETSC_FALSE;
  PetscFunctionReturn(0);
}

PETSC_STATIC_INLINE PetscErrorCode GetPointArray_Private(DM dm,PetscInt n,PetscInt *points,PetscInt *rn,const PetscInt **rpoints)
{
  PetscErrorCode ierr;
  PetscInt       *work;

  PetscFunctionBegin;
  if (rn) *rn = n;
  if (rpoints) {
    ierr     = DMGetWorkArray(dm,n,PETSC_INT,&work);CHKERRQ(ierr);
    ierr     = PetscMemcpy(work,points,n*sizeof(PetscInt));CHKERRQ(ierr);
    *rpoints = work;
  }
  PetscFunctionReturn(0);
}

/*@
   PetscCommSplitReductionBegin - Begin an asynchronous split-mode reduction

   Collective but not synchronizing

   Input Arguments:
   comm - communicator on which split reduction has been queued

   Level: advanced

   Note:
   Calling this function is optional when using split-mode reduction. On supporting hardware, calling this after all
   VecXxxBegin() allows the reduction to make asynchronous progress before the result is needed (in VecXxxEnd()).

.seealso: VecNormBegin(), VecNormEnd(), VecDotBegin(), VecDotEnd(), VecTDotBegin(), VecTDotEnd(), VecMDotBegin(), VecMDotEnd(), VecMTDotBegin(), VecMTDotEnd()
@*/
PetscErrorCode PetscCommSplitReductionBegin(MPI_Comm comm)
{
  PetscErrorCode      ierr;
  PetscSplitReduction *sr;

  PetscFunctionBegin;
  ierr = PetscSplitReductionGet(comm,&sr);CHKERRQ(ierr);
  if (sr->numopsend > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Cannot call this after VecxxxEnd() has been called");
  if (sr->async) {              /* Bad reuse, setup code copied from PetscSplitReductionApply(). */
    PetscInt       i,numops = sr->numopsbegin,*reducetype = sr->reducetype;
    PetscScalar    *lvalues = sr->lvalues,*gvalues = sr->gvalues;
    PetscInt       sum_flg = 0,max_flg = 0, min_flg = 0;
    MPI_Comm       comm = sr->comm;
    PetscMPIInt    size,cmul = sizeof(PetscScalar)/sizeof(PetscReal);;
    ierr = PetscLogEventBegin(VEC_ReduceBegin,0,0,0,0);CHKERRQ(ierr);
    ierr = MPI_Comm_size(sr->comm,&size);CHKERRQ(ierr);
    if (size == 1) {
      ierr = PetscMemcpy(gvalues,lvalues,numops*sizeof(PetscScalar));CHKERRQ(ierr);
    } else {
      /* determine if all reductions are sum, max, or min */
      for (i=0; i<numops; i++) {
        if      (reducetype[i] == REDUCE_MAX) max_flg = 1;
        else if (reducetype[i] == REDUCE_SUM) sum_flg = 1;
        else if (reducetype[i] == REDUCE_MIN) min_flg = 1;
        else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in PetscSplitReduction() data structure, probably memory corruption");
      }
      if (sum_flg + max_flg + min_flg > 1) {
        /*
         after all the entires in lvalues we store the reducetype flags to indicate
         to the reduction operations what are sums and what are max
         */
        for (i=0; i<numops; i++) lvalues[numops+i] = reducetype[i];

        ierr = MPIPetsc_Iallreduce(lvalues,gvalues,2*numops,MPIU_SCALAR,PetscSplitReduction_Op,comm,&sr->request);CHKERRQ(ierr);
      } else if (max_flg) {   /* Compute max of real and imag parts separately, presumably only the real part is used */
        ierr = MPIPetsc_Iallreduce((PetscReal*)lvalues,(PetscReal*)gvalues,cmul*numops,MPIU_REAL,MPIU_MAX,comm,&sr->request);CHKERRQ(ierr);
      } else if (min_flg) {
        ierr = MPIPetsc_Iallreduce((PetscReal*)lvalues,(PetscReal*)gvalues,cmul*numops,MPIU_REAL,MPIU_MIN,comm,&sr->request);CHKERRQ(ierr);
      } else {
        ierr = MPIPetsc_Iallreduce(lvalues,gvalues,numops,MPIU_SCALAR,MPIU_SUM,comm,&sr->request);CHKERRQ(ierr);
      }
    }
    sr->state     = STATE_PENDING;
    sr->numopsend = 0;
    ierr = PetscLogEventEnd(VEC_ReduceBegin,0,0,0,0);CHKERRQ(ierr);
  } else {
    ierr = PetscSplitReductionApply(sr);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

PetscErrorCode MatHeaderMerge(Mat A,Mat C)
{
  PetscErrorCode ierr;
  PetscInt       refct;
  PetscOps       *Abops;
  MatOps         Aops;
  char           *mtype,*mname;
  void           *spptr;

  PetscFunctionBegin;
  /* save the parts of A we need */
  Abops = ((PetscObject)A)->bops;
  Aops  = A->ops;
  refct = ((PetscObject)A)->refct;
  mtype = ((PetscObject)A)->type_name;
  mname = ((PetscObject)A)->name;
  spptr = A->spptr;

  /* zero these so the destroy below does not free them */
  ((PetscObject)A)->type_name = 0;
  ((PetscObject)A)->name      = 0;

  /* free all the interior data structures from mat */
  ierr = (*A->ops->destroy)(A);CHKERRQ(ierr);

  ierr = PetscFree(C->spptr);CHKERRQ(ierr);

  ierr = PetscLayoutDestroy(&A->rmap);CHKERRQ(ierr);
  ierr = PetscLayoutDestroy(&A->cmap);CHKERRQ(ierr);
  ierr = PetscFunctionListDestroy(&((PetscObject)A)->qlist);CHKERRQ(ierr);
  ierr = PetscObjectListDestroy(&((PetscObject)A)->olist);CHKERRQ(ierr);

  /* copy C over to A */
  ierr = PetscMemcpy(A,C,sizeof(struct _p_Mat));CHKERRQ(ierr);

  /* return the parts of A we saved */
  ((PetscObject)A)->bops      = Abops;
  A->ops                      = Aops;
  ((PetscObject)A)->refct     = refct;
  ((PetscObject)A)->type_name = mtype;
  ((PetscObject)A)->name      = mname;
  A->spptr                    = spptr;

  /* since these two are copied into A we do not want them destroyed in C */
  ((PetscObject)C)->qlist = 0;
  ((PetscObject)C)->olist = 0;

  ierr = PetscHeaderDestroy(&C);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*@
   PetscMergeIntArray -     Merges two SORTED integer arrays, removes duplicate elements.

   Not Collective

   Input Parameters:
+  an  - number of values in the first array
.  aI  - first sorted array of integers
.  bn  - number of values in the second array
-  bI  - second array of integers

   Output Parameters:
+  n   - number of values in the merged array
-  L   - merged sorted array, this is allocated if an array is not provided

   Level: intermediate

   Concepts: merging^arrays

.seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
@*/
PetscErrorCode  PetscMergeIntArray(PetscInt an,const PetscInt aI[], PetscInt bn, const PetscInt bI[],  PetscInt *n, PetscInt **L)
{
  PetscErrorCode ierr;
  PetscInt       *L_ = *L, ak, bk, k;

  if (!L_) {
    ierr = PetscMalloc1(an+bn, L);CHKERRQ(ierr);
    L_   = *L;
  }
  k = ak = bk = 0;
  while (ak < an && bk < bn) {
    if (aI[ak] == bI[bk]) {
      L_[k] = aI[ak];
      ++ak;
      ++bk;
      ++k;
    } else if (aI[ak] < bI[bk]) {
      L_[k] = aI[ak];
      ++ak;
      ++k;
    } else {
      L_[k] = bI[bk];
      ++bk;
      ++k;
    }
  }
  if (ak < an) {
    ierr = PetscMemcpy(L_+k,aI+ak,(an-ak)*sizeof(PetscInt));CHKERRQ(ierr);
    k   += (an-ak);
  }
  if (bk < bn) {
    ierr = PetscMemcpy(L_+k,bI+bk,(bn-bk)*sizeof(PetscInt));CHKERRQ(ierr);
    k   += (bn-bk);
  }
  *n = k;
  PetscFunctionReturn(0);
}

/*@C
   TSRKRegister - register an RK scheme by providing the entries in the Butcher tableau and optionally embedded approximations and interpolation

   Not Collective, but the same schemes should be registered on all processes on which they will be used

   Input Parameters:
+  name - identifier for method
.  order - approximation order of method
.  s - number of stages, this is the dimension of the matrices below
.  A - stage coefficients (dimension s*s, row-major)
.  b - step completion table (dimension s; NULL to use last row of A)
.  c - abscissa (dimension s; NULL to use row sums of A)
.  bembed - completion table for embedded method (dimension s; NULL if not available)
.  pinterp - Order of the interpolation scheme, equal to the number of columns of binterp
-  binterp - Coefficients of the interpolation formula (dimension s*pinterp; NULL to reuse binterpt)

   Notes:
   Several RK methods are provided, this function is only needed to create new methods.

   Level: advanced

.keywords: TS, register

.seealso: TSRK
@*/
PetscErrorCode TSRKRegister(TSRKType name,PetscInt order,PetscInt s,
                                 const PetscReal A[],const PetscReal b[],const PetscReal c[],
                                 const PetscReal bembed[],
                                 PetscInt pinterp,const PetscReal binterp[])
{
  PetscErrorCode  ierr;
  RKTableauLink   link;
  RKTableau       t;
  PetscInt        i,j;

  PetscFunctionBegin;
  ierr     = PetscMalloc(sizeof(*link),&link);CHKERRQ(ierr);
  ierr     = PetscMemzero(link,sizeof(*link));CHKERRQ(ierr);
  t        = &link->tab;
  ierr     = PetscStrallocpy(name,&t->name);CHKERRQ(ierr);
  t->order = order;
  t->s     = s;
  ierr     = PetscMalloc3(s*s,&t->A,s,&t->b,s,&t->c);CHKERRQ(ierr);
  ierr     = PetscMemcpy(t->A,A,s*s*sizeof(A[0]));CHKERRQ(ierr);
  if (b)  { ierr = PetscMemcpy(t->b,b,s*sizeof(b[0]));CHKERRQ(ierr); }
  else for (i=0; i<s; i++) t->b[i] = A[(s-1)*s+i];
  if (c)  { ierr = PetscMemcpy(t->c,c,s*sizeof(c[0]));CHKERRQ(ierr); }
  else for (i=0; i<s; i++) for (j=0,t->c[i]=0; j<s; j++) t->c[i] += A[i*s+j];
  t->FSAL = PETSC_TRUE;
  for (i=0; i<s; i++) if (t->A[(s-1)*s+i] != t->b[i]) t->FSAL = PETSC_FALSE;
  if (bembed) {
    ierr = PetscMalloc1(s,&t->bembed);CHKERRQ(ierr);
    ierr = PetscMemcpy(t->bembed,bembed,s*sizeof(bembed[0]));CHKERRQ(ierr);
  }

  t->pinterp     = pinterp;
  ierr           = PetscMalloc1(s*pinterp,&t->binterp);CHKERRQ(ierr);
  ierr           = PetscMemcpy(t->binterp,binterp,s*pinterp*sizeof(binterp[0]));CHKERRQ(ierr);
  link->next     = RKTableauList;
  RKTableauList = link;
  PetscFunctionReturn(0);
}

/*@C
  ADDAHCiterStartup - performs the first check for an iteration through a hypercube
  lc, uc, idx all have to be valid arrays of size dim
  This function sets idx to lc and then checks, whether the lower corner (lc) is less
  than thre upper corner (uc). If lc "<=" uc in all coordinates, it returns PETSC_TRUE,
  and PETSC_FALSE otherwise.

  Input Parameters:
+ dim - the number of dimension
. lc - the "lower" corner
- uc - the "upper" corner

  Output Parameters:
. idx - the index that this function increases

  Developer Notes: This code is crap! You cannot return a value and NO ERROR code in PETSc!

  Level: developer
@*/
PetscBool  ADDAHCiterStartup(const PetscInt dim, const PetscInt *const lc, const PetscInt *const uc, PetscInt *const idx)
{
  PetscErrorCode ierr;
  PetscInt       i;

  ierr = PetscMemcpy(idx, lc, sizeof(PetscInt)*dim);
  if (ierr) {
    PetscError(PETSC_COMM_SELF,__LINE__,__FUNCT__,__FILE__,__SDIR__,ierr,PETSC_ERROR_REPEAT," ");
    return PETSC_FALSE;
  }
  for (i=0; i<dim; i++) {
    if (lc[i] > uc[i]) return PETSC_FALSE;
  }
  return PETSC_TRUE;
}

/*@C
    DMSlicedSetGhosts - Sets the global indices of other processes elements that will
      be ghosts on this process

    Not Collective

    Input Parameters:
+    slice - the DM object
.    bs - block size
.    nlocal - number of local (owned, non-ghost) blocks
.    Nghosts - number of ghost blocks on this process
-    ghosts - global indices of each ghost block

    Level: advanced

.seealso DMDestroy(), DMCreateGlobalVector()

@*/
PetscErrorCode  DMSlicedSetGhosts(DM dm,PetscInt bs,PetscInt nlocal,PetscInt Nghosts,const PetscInt ghosts[])
{
  PetscErrorCode ierr;
  DM_Sliced      *slice = (DM_Sliced*)dm->data;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
  ierr = PetscFree(slice->ghosts);CHKERRQ(ierr);
  ierr = PetscMalloc(Nghosts*sizeof(PetscInt),&slice->ghosts);CHKERRQ(ierr);
  ierr = PetscMemcpy(slice->ghosts,ghosts,Nghosts*sizeof(PetscInt));CHKERRQ(ierr);
  slice->bs      = bs;
  slice->n       = nlocal;
  slice->Nghosts = Nghosts;
  PetscFunctionReturn(0);
}

PetscErrorCode PetscDrawSetColormap_X(PetscDraw_X *XiWin,Colormap colormap)
{
  PetscErrorCode ierr;

  PetscFunctionBegin;
  if (XiWin->depth < 8) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP_SYS,"PETSc Graphics require monitors with at least 8 bit color (256 colors)");
  if (!gColormap) {
    ierr = PetscDrawSetUpColormap_X(XiWin->disp,XiWin->screen,XiWin->vis,colormap);CHKERRQ(ierr);
  }
  XiWin->cmap       = gColormap;
  ierr              = PetscMemcpy(XiWin->cmapping,gCmapping,256*sizeof(PetscDrawXiPixVal));CHKERRQ(ierr);
  XiWin->background = XiWin->cmapping[PETSC_DRAW_WHITE];
  XiWin->foreground = XiWin->cmapping[PETSC_DRAW_BLACK];
  PetscFunctionReturn(0);
}

EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "PetscRandomCreate_Rand48"
PetscErrorCode  PetscRandomCreate_Rand48(PetscRandom r)
{
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = PetscMemcpy(r->ops,&PetscRandomOps_Values,sizeof(PetscRandomOps_Values));CHKERRQ(ierr);
  /* r->bops->publish   = PetscRandomPublish; */
  /*  r->petscnative     = PETSC_TRUE;  */

  ierr = PetscObjectChangeTypeName((PetscObject)r,PETSCRAND48);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode  ISBlockSetIndices_Block(IS is,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
{
  PetscErrorCode ierr;
  PetscInt       i,min,max;
  IS_Block       *sub   = (IS_Block*)is->data;
  PetscBool      sorted = PETSC_TRUE;

  PetscFunctionBegin;
  ierr = PetscFree(sub->idx);CHKERRQ(ierr);
  sub->n = n;
  ierr = MPI_Allreduce(&n,&sub->N,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)is));CHKERRQ(ierr);
  for (i=1; i<n; i++) {
    if (idx[i] < idx[i-1]) {sorted = PETSC_FALSE; break;}
  }
  if (n) min = max = idx[0];
  else   min = max = 0;
  for (i=1; i<n; i++) {
    if (idx[i] < min) min = idx[i];
    if (idx[i] > max) max = idx[i];
  }
  if (mode == PETSC_COPY_VALUES) {
    ierr = PetscMalloc(n*sizeof(PetscInt),&sub->idx);CHKERRQ(ierr);
    ierr = PetscLogObjectMemory(is,n*sizeof(PetscInt));CHKERRQ(ierr);
    ierr = PetscMemcpy(sub->idx,idx,n*sizeof(PetscInt));CHKERRQ(ierr);
  } else if (mode == PETSC_OWN_POINTER) sub->idx = (PetscInt*) idx;
  else SETERRQ(PetscObjectComm((PetscObject)is),PETSC_ERR_SUP,"Only supports PETSC_COPY_VALUES and PETSC_OWN_POINTER");

  sub->sorted = sorted;
  is->bs      = bs;
  is->min     = bs*min;
  is->max     = bs*max+bs-1;
  is->data    = (void*)sub;
  ierr = PetscMemcpy(is->ops,&myops,sizeof(myops));CHKERRQ(ierr);
  is->isperm  = PETSC_FALSE;
  PetscFunctionReturn(0);
}

/*@C
   ISGetNonlocalIndices - Retrieve an array of indices from remote processors
                       in this communicator.

   Collective on IS

   Input Parameter:
.  is - the index set

   Output Parameter:
.  indices - indices with rank 0 indices first, and so on,  omitting
             the current rank.  Total number of indices is the difference
             total and local, obtained with ISGetSize() and ISGetLocalSize(),
             respectively.

   Level: intermediate

   Notes: restore the indices using ISRestoreNonlocalIndices().
          The same scalability considerations as those for ISGetTotalIndices
          apply here.

   Concepts: index sets^getting nonlocal indices
.seealso: ISGetTotalIndices(), ISRestoreNonlocalIndices(), ISGetSize(), ISGetLocalSize().
@*/
PetscErrorCode  ISGetNonlocalIndices(IS is, const PetscInt *indices[])
{
  PetscErrorCode ierr;
  PetscMPIInt    size;
  PetscInt       n, N;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(is,IS_CLASSID,1);
  PetscValidPointer(indices,2);
  ierr = MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);CHKERRQ(ierr);
  if (size == 1) *indices = NULL;
  else {
    if (!is->total) {
      ierr = ISGatherTotal_Private(is);CHKERRQ(ierr);
    }
    ierr = ISGetLocalSize(is,&n);CHKERRQ(ierr);
    ierr = ISGetSize(is,&N);CHKERRQ(ierr);
    ierr = PetscMalloc(sizeof(PetscInt)*(N-n), &(is->nonlocal));CHKERRQ(ierr);
    ierr = PetscMemcpy(is->nonlocal, is->total, sizeof(PetscInt)*is->local_offset);CHKERRQ(ierr);
    ierr = PetscMemcpy(is->nonlocal+is->local_offset, is->total+is->local_offset+n, sizeof(PetscInt)*(N - is->local_offset - n));CHKERRQ(ierr);
    *indices = is->nonlocal;
  }
  PetscFunctionReturn(0);
}

PetscErrorCode SNESVIRedundancyCheck_Matlab(SNES snes,IS is_act,IS* is_redact,void* ctx)
{
  PetscErrorCode      ierr;
  SNESMatlabContext   *sctx = (SNESMatlabContext*)ctx;
  int                 nlhs = 1, nrhs = 5;
  mxArray             *plhs[1], *prhs[5];
  long long int       l1 = 0, l2 = 0, ls = 0;
  PetscInt            *indices=PETSC_NULL;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
  PetscValidHeaderSpecific(is_act,IS_CLASSID,2);
  PetscValidPointer(is_redact,3);
  PetscCheckSameComm(snes,1,is_act,2);

  /* Create IS for reduced active set of size 0, its size and indices will
   bet set by the Matlab function */
  ierr = ISCreateGeneral(((PetscObject)snes)->comm,0,indices,PETSC_OWN_POINTER,is_redact);CHKERRQ(ierr);
  /* call Matlab function in ctx */
  ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr);
  ierr = PetscMemcpy(&l1,&is_act,sizeof(is_act));CHKERRQ(ierr);
  ierr = PetscMemcpy(&l2,is_redact,sizeof(is_act));CHKERRQ(ierr);
  prhs[0] = mxCreateDoubleScalar((double)ls);
  prhs[1] = mxCreateDoubleScalar((double)l1);
  prhs[2] = mxCreateDoubleScalar((double)l2);
  prhs[3] = mxCreateString(sctx->funcname);
  prhs[4] = sctx->ctx;
  ierr    = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESVIRedundancyCheckInternal");CHKERRQ(ierr);
  ierr    = mxGetScalar(plhs[0]);CHKERRQ(ierr);
  mxDestroyArray(prhs[0]);
  mxDestroyArray(prhs[1]);
  mxDestroyArray(prhs[2]);
  mxDestroyArray(prhs[3]);
  mxDestroyArray(plhs[0]);
  PetscFunctionReturn(0);
}

/*@C
   PetscFortranCallbackRegister - register a type+subtype callback

   Not Collective

   Input Arguments:
+  classid - ID of class on which to register callback
-  subtype - subtype string, or NULL for class ids

   Output Arguments:
.  id - callback id

   Level: developer

.seealso: PetscFortranCallbackGetSizes()
@*/
PetscErrorCode PetscFortranCallbackRegister(PetscClassId classid,const char *subtype,PetscFortranCallbackId *id)
{
  PetscErrorCode      ierr;
  FortranCallbackBase *base;
  FortranCallbackLink link;

  PetscFunctionBegin;
  *id = 0;
  if (classid < PETSC_SMALLEST_CLASSID || PETSC_LARGEST_CLASSID <= classid) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"ClassId %D corrupt",classid);
  if (classid >= _maxclassid) {
    PetscClassId        newmax = PETSC_SMALLEST_CLASSID + 2*(PETSC_LARGEST_CLASSID-PETSC_SMALLEST_CLASSID);
    FortranCallbackBase *newbase;
    if (!_classbase) {
      ierr = PetscRegisterFinalize(PetscFortranCallbackFinalize);CHKERRQ(ierr);
    }
    ierr = PetscMalloc((newmax-PETSC_SMALLEST_CLASSID)*sizeof(_classbase[0]),&newbase);CHKERRQ(ierr);
    ierr = PetscMemzero(newbase,(newmax-PETSC_SMALLEST_CLASSID)*sizeof(_classbase[0]));CHKERRQ(ierr);
    ierr = PetscMemcpy(newbase,_classbase,(_maxclassid-PETSC_SMALLEST_CLASSID)*sizeof(_classbase[0]));CHKERRQ(ierr);
    ierr = PetscFree(_classbase);CHKERRQ(ierr);

    _classbase = newbase;
    _maxclassid = newmax;
  }
  base = &_classbase[classid-PETSC_SMALLEST_CLASSID];
  if (!subtype) *id = PETSC_SMALLEST_FORTRAN_CALLBACK + base->basecount++;
  else {
    for (link=base->subtypes; link; link=link->next) { /* look for either both NULL or matching values (implies both non-NULL) */
      PetscBool match;
      ierr = PetscStrcmp(subtype,link->type_name,&match);CHKERRQ(ierr);
      if (match) { /* base type or matching subtype */
        goto found;
      }
    }
    /* Not found. Create node and prepend to class' subtype list */
    ierr = PetscMalloc(sizeof(*link),&link);CHKERRQ(ierr);
    ierr = PetscStrallocpy(subtype,&link->type_name);CHKERRQ(ierr);

    link->max      = PETSC_SMALLEST_FORTRAN_CALLBACK;
    link->next     = base->subtypes;
    base->subtypes = link;

found:
    *id = link->max++;

    base->maxsubtypecount = PetscMax(base->maxsubtypecount,link->max-PETSC_SMALLEST_FORTRAN_CALLBACK);
  }
  PetscFunctionReturn(0);
}

PetscErrorCode BVCopy_Mat(BV V,BV W)
{
  PetscErrorCode ierr;
  BV_MAT         *v = (BV_MAT*)V->data,*w = (BV_MAT*)W->data;
  PetscScalar    *pv,*pw,*pvc,*pwc;

  PetscFunctionBegin;
  ierr = MatDenseGetArray(v->A,&pv);CHKERRQ(ierr);
  ierr = MatDenseGetArray(w->A,&pw);CHKERRQ(ierr);
  pvc = pv+(V->nc+V->l)*V->n;
  pwc = pw+(W->nc+W->l)*W->n;
  ierr = PetscMemcpy(pwc,pvc,(V->k-V->l)*V->n*sizeof(PetscScalar));CHKERRQ(ierr);
  ierr = MatDenseRestoreArray(v->A,&pv);CHKERRQ(ierr);
  ierr = MatDenseRestoreArray(w->A,&pw);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode PetscViewerGetSingleton_Binary(PetscViewer viewer,PetscViewer *outviewer)
{
  int                rank;
  PetscErrorCode     ierr;
  PetscViewer_Binary *vbinary = (PetscViewer_Binary*)viewer->data,*obinary;

  PetscFunctionBegin;
  ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)viewer),&rank);CHKERRQ(ierr);
  if (!rank) {
    ierr    = PetscViewerCreate(PETSC_COMM_SELF,outviewer);CHKERRQ(ierr);
    ierr    = PetscViewerSetType(*outviewer,PETSCVIEWERBINARY);CHKERRQ(ierr);
    obinary = (PetscViewer_Binary*)(*outviewer)->data;
    ierr    = PetscMemcpy(obinary,vbinary,sizeof(PetscViewer_Binary));CHKERRQ(ierr);
  } else *outviewer = 0;
  PetscFunctionReturn(0);
}

PETSC_EXTERN PetscMPIInt PetscDataRep_read_conv_fn(void *userbuf, MPI_Datatype datatype,PetscMPIInt count,void *filebuf, MPI_Offset position,void *extra_state)
{
  PetscDataType pdtype;
  PetscMPIInt   ierr;
  size_t        dsize;

  ierr = PetscMPIDataTypeToPetscDataType(datatype,&pdtype);CHKERRQ(ierr);
  ierr = PetscDataTypeGetSize(pdtype,&dsize);CHKERRQ(ierr);

  /* offset is given in units of MPI_Datatype */
  userbuf = ((char*)userbuf) + dsize*position;

  ierr = PetscMemcpy(userbuf,filebuf,count*dsize);CHKERRQ(ierr);
  ierr = PetscByteSwap(userbuf,pdtype,count);CHKERRQ(ierr);
  return ierr;
}