PetscErrorCode MatMPIAIJDiagonalScaleLocalSetUp(Mat inA,Vec scale)
{
  Mat_MPIAIJ     *ina = (Mat_MPIAIJ*) inA->data; /*access private part of matrix */
  PetscErrorCode ierr;
  PetscInt       i,n,nt,cstart,cend,no,*garray = ina->garray,*lindices;
  PetscInt       *r_rmapd,*r_rmapo;

  PetscFunctionBegin;
  ierr = MatGetOwnershipRange(inA,&cstart,&cend);CHKERRQ(ierr);
  ierr = MatGetSize(ina->A,NULL,&n);CHKERRQ(ierr);
  ierr = PetscCalloc1(inA->rmap->mapping->n+1,&r_rmapd);CHKERRQ(ierr);
  nt   = 0;
  for (i=0; i<inA->rmap->mapping->n; i++) {
    if (inA->rmap->mapping->indices[i] >= cstart && inA->rmap->mapping->indices[i] < cend) {
      nt++;
      r_rmapd[i] = inA->rmap->mapping->indices[i] + 1;
    }
  }
  if (nt != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Hmm nt %D n %D",nt,n);
  ierr = PetscMalloc1(n+1,&auglyrmapd);CHKERRQ(ierr);
  for (i=0; i<inA->rmap->mapping->n; i++) {
    if (r_rmapd[i]) {
      auglyrmapd[(r_rmapd[i]-1)-cstart] = i;
    }
  }
  ierr = PetscFree(r_rmapd);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF,n,&auglydd);CHKERRQ(ierr);

  ierr = PetscCalloc1(inA->cmap->N+1,&lindices);CHKERRQ(ierr);
  for (i=0; i<ina->B->cmap->n; i++) {
    lindices[garray[i]] = i+1;
  }
  no   = inA->rmap->mapping->n - nt;
  ierr = PetscCalloc1(inA->rmap->mapping->n+1,&r_rmapo);CHKERRQ(ierr);
  nt   = 0;
  for (i=0; i<inA->rmap->mapping->n; i++) {
    if (lindices[inA->rmap->mapping->indices[i]]) {
      nt++;
      r_rmapo[i] = lindices[inA->rmap->mapping->indices[i]];
    }
  }
  if (nt > no) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Hmm nt %D no %D",nt,n);
  ierr = PetscFree(lindices);CHKERRQ(ierr);
  ierr = PetscMalloc1(nt+1,&auglyrmapo);CHKERRQ(ierr);
  for (i=0; i<inA->rmap->mapping->n; i++) {
    if (r_rmapo[i]) {
      auglyrmapo[(r_rmapo[i]-1)] = i;
    }
  }
  ierr = PetscFree(r_rmapo);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF,nt,&auglyoo);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode MatPartitioningHierarchical_ReassembleFineparts(Mat adj, IS fineparts, ISLocalToGlobalMapping mapping, IS *sfineparts)
{
  PetscInt            *local_indices, *global_indices,*owners,*sfineparts_indices,localsize,i;
  const PetscInt      *ranges,*fineparts_indices;
  PetscMPIInt         rank;
  MPI_Comm            comm;
  PetscLayout         rmap;
  PetscSFNode        *remote;
  PetscSF             sf;
  PetscErrorCode      ierr;

  PetscFunctionBegin;
  /*get communicator */
  ierr = PetscObjectGetComm((PetscObject)adj,&comm);CHKERRQ(ierr);
  ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
  ierr = MatGetLayouts(adj,&rmap,PETSC_NULL);CHKERRQ(ierr);
  ierr = ISGetLocalSize(fineparts,&localsize);CHKERRQ(ierr);
  ierr = PetscCalloc2(localsize,&global_indices,localsize,&local_indices);CHKERRQ(ierr);
  for(i=0; i<localsize; i++){
	local_indices[i] = i;
  }
  /*global indices */
  ierr = ISLocalToGlobalMappingApply(mapping,localsize,local_indices,global_indices);CHKERRQ(ierr);
  ierr = PetscCalloc1(localsize,&owners);CHKERRQ(ierr);
  /*find owners for global indices */
  for(i=0; i<localsize; i++){
	ierr = PetscLayoutFindOwner(rmap,global_indices[i],&owners[i]);CHKERRQ(ierr);
  }
  /*ranges */
  ierr = PetscLayoutGetRanges(rmap,&ranges);CHKERRQ(ierr);
  ierr = PetscCalloc1(ranges[rank+1]-ranges[rank],&sfineparts_indices);CHKERRQ(ierr);
  ierr = ISGetIndices(fineparts,&fineparts_indices);CHKERRQ(ierr);
  /*create a SF to exchange data */
  ierr = PetscSFCreate(comm,&sf);CHKERRQ(ierr);
  ierr = PetscCalloc1(localsize,&remote);CHKERRQ(ierr);
  for(i=0; i<localsize; i++){
	remote[i].rank  = owners[i];
	remote[i].index = global_indices[i]-ranges[owners[i]];
  }
  ierr = PetscSFSetType(sf,PETSCSFBASIC);CHKERRQ(ierr);
  /*not sure how to add prefix to sf*/
  ierr = PetscSFSetFromOptions(sf);CHKERRQ(ierr);
  ierr = PetscSFSetGraph(sf,localsize,localsize,PETSC_NULL,PETSC_OWN_POINTER,remote,PETSC_OWN_POINTER);CHKERRQ(ierr);
  ierr = PetscSFReduceBegin(sf,MPIU_INT,fineparts_indices,sfineparts_indices,MPIU_REPLACE);CHKERRQ(ierr);
  ierr = PetscSFReduceEnd(sf,MPIU_INT,fineparts_indices,sfineparts_indices,MPIU_REPLACE);CHKERRQ(ierr);
  ierr = PetscSFDestroy(&sf);CHKERRQ(ierr);
  ierr = ISRestoreIndices(fineparts,&fineparts_indices);CHKERRQ(ierr);
  /* comm self */
  ierr = ISCreateGeneral(comm,ranges[rank+1]-ranges[rank],sfineparts_indices,PETSC_OWN_POINTER,sfineparts);CHKERRQ(ierr);
  ierr = PetscFree2(global_indices,local_indices);CHKERRQ(ierr);
  ierr = PetscFree(owners);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

static PetscErrorCode CreatePoints_Grid(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx)
{
  DM             da;
  DMDALocalInfo  info;
  PetscInt       N = 3, n = 0, spaceDim, i, j, k, *ind, d;
  PetscReal      *h;
  PetscMPIInt    rank;
  PetscErrorCode ierr;

  ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRQ(ierr);
  ierr = DMGetCoordinateDim(dm, &spaceDim);CHKERRQ(ierr);
  ierr = PetscCalloc1(spaceDim,&ind);CHKERRQ(ierr);
  ierr = PetscCalloc1(spaceDim,&h);CHKERRQ(ierr);
  h[0] = 1.0/(N-1); h[1] = 1.0/(N-1); h[2] = 1.0/(N-1);
  ierr = DMDACreate(PetscObjectComm((PetscObject) dm), &da);CHKERRQ(ierr);
  ierr = DMSetDimension(da, ctx->dim);CHKERRQ(ierr);
  ierr = DMDASetSizes(da, N, N, N);CHKERRQ(ierr);
  ierr = DMDASetDof(da, 1);CHKERRQ(ierr);
  ierr = DMDASetStencilWidth(da, 1);CHKERRQ(ierr);
  ierr = DMSetUp(da);CHKERRQ(ierr);
  ierr = DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);CHKERRQ(ierr);
  ierr = DMDAGetLocalInfo(da, &info);CHKERRQ(ierr);
  *Np  = info.xm * info.ym * info.zm;
  ierr = PetscCalloc1(*Np * spaceDim, pcoords);CHKERRQ(ierr);
  for (k = info.zs; k < info.zs + info.zm; ++k) {
    ind[2] = k;
    for (j = info.ys; j < info.ys + info.ym; ++j) {
      ind[1] = j;
      for (i = info.xs; i < info.xs + info.xm; ++i, ++n) {
        ind[0] = i;

        for (d = 0; d < spaceDim; ++d) (*pcoords)[n*spaceDim+d] = ind[d]*h[d];
        ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %D (", rank, n);CHKERRQ(ierr);
        for (d = 0; d < spaceDim; ++d) {
          ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%g", (double)(*pcoords)[n*spaceDim+d]);CHKERRQ(ierr);
          if (d < spaceDim-1) {ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, ", ");CHKERRQ(ierr);}
        }
        ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, ")\n");CHKERRQ(ierr);
      }
    }
  }
  ierr = DMDestroy(&da);CHKERRQ(ierr);
  ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL);CHKERRQ(ierr);
  ierr = PetscFree(ind);CHKERRQ(ierr);
  ierr = PetscFree(h);CHKERRQ(ierr);
  *pointsAllProcs = PETSC_FALSE;
  PetscFunctionReturn(0);
}

/*@C
  DMPlexCreateLabel - Create a label of the given name if it does not already exist

  Not Collective

  Input Parameters:
+ dm   - The DMPlex object
- name - The label name

  Level: intermediate

.keywords: mesh
.seealso: DMLabelCreate(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
@*/
PetscErrorCode DMPlexCreateLabel(DM dm, const char name[])
{
  DM_Plex        *mesh = (DM_Plex*) dm->data;
  PlexLabel      next  = mesh->labels;
  PetscBool      flg   = PETSC_FALSE;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
  PetscValidCharPointer(name, 2);
  while (next) {
    ierr = PetscStrcmp(name, next->label->name, &flg);CHKERRQ(ierr);
    if (flg) break;
    next = next->next;
  }
  if (!flg) {
    PlexLabel tmpLabel;

    ierr = PetscCalloc1(1, &tmpLabel);CHKERRQ(ierr);
    ierr = DMLabelCreate(name, &tmpLabel->label);CHKERRQ(ierr);
    tmpLabel->output = PETSC_TRUE;
    tmpLabel->next   = mesh->labels;
    mesh->labels     = tmpLabel;
  }
  PetscFunctionReturn(0);
}

/*
    PetscTableAddExpand - called PetscTableAdd() if more space needed

*/
PetscErrorCode  PetscTableAddExpand(PetscTable ta,PetscInt key,PetscInt data,InsertMode imode)
{
  PetscErrorCode ierr;
  PetscInt       ii      = 0;
  const PetscInt tsize   = ta->tablesize,tcount = ta->count;
  PetscInt       *oldtab = ta->table,*oldkt = ta->keytable,newk,ndata;

  PetscFunctionBegin;
  if (ta->tablesize == PETSC_MAX_INT/4) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_COR,"ta->tablesize < 0");
  ta->tablesize = 2*tsize;
  if (ta->tablesize <= tsize) ta->tablesize = PETSC_MAX_INT/4;

  ierr = PetscMalloc1(ta->tablesize,&ta->table);CHKERRQ(ierr);
  ierr = PetscCalloc1(ta->tablesize,&ta->keytable);CHKERRQ(ierr);

  ta->count = 0;
  ta->head  = 0;

  ierr = PetscTableAdd(ta,key,data,INSERT_VALUES);CHKERRQ(ierr);
  /* rehash */
  for (ii = 0; ii < tsize; ii++) {
    newk = oldkt[ii];
    if (newk) {
      ndata = oldtab[ii];
      ierr  = PetscTableAdd(ta,newk,ndata,imode);CHKERRQ(ierr);
    }
  }
  if (ta->count != tcount + 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_COR,"corrupted ta->count");

  ierr = PetscFree(oldtab);CHKERRQ(ierr);
  ierr = PetscFree(oldkt);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode MatPartitioningHierarchical_DetermineDestination(MatPartitioning part, IS partitioning, PetscInt pstart, PetscInt pend, IS *destination)
{
  MPI_Comm            comm;
  PetscMPIInt         rank,size,target;
  PetscInt            plocalsize,*dest_indices,i;
  const PetscInt     *part_indices;
  PetscErrorCode      ierr;

  PetscFunctionBegin;
  /*communicator*/
  ierr = PetscObjectGetComm((PetscObject)part,&comm);CHKERRQ(ierr);
  ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
  ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
  if((pend-pstart)>size) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"range [%D, %D] should be smaller than or equal to size %D",pstart,pend,size);CHKERRQ(ierr);
  if(pstart>pend) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP," pstart %D should be smaller than pend %D",pstart,pend);CHKERRQ(ierr);
  /*local size*/
  ierr = ISGetLocalSize(partitioning,&plocalsize);CHKERRQ(ierr);
  ierr = PetscCalloc1(plocalsize,&dest_indices);CHKERRQ(ierr);
  ierr = ISGetIndices(partitioning,&part_indices);CHKERRQ(ierr);
  for(i=0; i<plocalsize; i++){
	/*compute target */
    target = part_indices[i]-pstart;
    /*mark out of range entity as -1*/
    if(part_indices[i]<pstart || part_indices[i]>pend) target = -1;
	dest_indices[i] = target;
  }
  /*return destination back*/
  ierr = ISCreateGeneral(comm,plocalsize,dest_indices,PETSC_OWN_POINTER,destination);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*@C
   PetscObjectCopyFortranFunctionPointers - Copy function pointers to another object

   Logically Collective on PetscObject

   Input Parameter:
+  src - source object
-  dest - destination object

   Level: developer

   Note:
   Both objects must have the same class.
@*/
PetscErrorCode PetscObjectCopyFortranFunctionPointers(PetscObject src,PetscObject dest)
{
  PetscErrorCode ierr;
  PetscInt       cbtype,numcb[PETSC_FORTRAN_CALLBACK_MAXTYPE];

  PetscFunctionBegin;
  PetscValidHeader(src,1);
  PetscValidHeader(dest,2);
  if (src->classid != dest->classid) SETERRQ(src->comm,PETSC_ERR_ARG_INCOMP,"Objects must be of the same class");

  ierr = PetscFree(dest->fortran_func_pointers);CHKERRQ(ierr);
  ierr = PetscMalloc(src->num_fortran_func_pointers*sizeof(void(*)(void)),&dest->fortran_func_pointers);CHKERRQ(ierr);
  ierr = PetscMemcpy(dest->fortran_func_pointers,src->fortran_func_pointers,src->num_fortran_func_pointers*sizeof(void(*)(void)));CHKERRQ(ierr);

  dest->num_fortran_func_pointers = src->num_fortran_func_pointers;

  ierr = PetscFortranCallbackGetSizes(src->classid,&numcb[PETSC_FORTRAN_CALLBACK_CLASS],&numcb[PETSC_FORTRAN_CALLBACK_SUBTYPE]);CHKERRQ(ierr);
  for (cbtype=PETSC_FORTRAN_CALLBACK_CLASS; cbtype<PETSC_FORTRAN_CALLBACK_MAXTYPE; cbtype++) {
    ierr = PetscFree(dest->fortrancallback[cbtype]);CHKERRQ(ierr);
    ierr = PetscCalloc1(numcb[cbtype],&dest->fortrancallback[cbtype]);CHKERRQ(ierr);
    ierr = PetscMemcpy(dest->fortrancallback[cbtype],src->fortrancallback[cbtype],src->num_fortrancallback[cbtype]*sizeof(PetscFortranCallback));CHKERRQ(ierr);
    dest->num_fortrancallback[cbtype] = src->num_fortrancallback[cbtype];
  }
  PetscFunctionReturn(0);
}

PetscErrorCode DMLabelDistribute_Internal(DMLabel label, PetscSF sf, PetscSection *leafSection, PetscInt **leafStrata)
{
  MPI_Comm       comm;
  PetscInt       s, l, nroots, nleaves, dof, offset, size;
  PetscInt      *remoteOffsets, *rootStrata, *rootIdx;
  PetscSection   rootSection;
  PetscSF        labelSF;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  if (label) {ierr = DMLabelMakeAllValid_Private(label);CHKERRQ(ierr);}
  ierr = PetscObjectGetComm((PetscObject)sf, &comm);CHKERRQ(ierr);
  /* Build a section of stratum values per point, generate the according SF
     and distribute point-wise stratum values to leaves. */
  ierr = PetscSFGetGraph(sf, &nroots, &nleaves, NULL, NULL);CHKERRQ(ierr);
  ierr = PetscSectionCreate(comm, &rootSection);CHKERRQ(ierr);
  ierr = PetscSectionSetChart(rootSection, 0, nroots);CHKERRQ(ierr);
  if (label) {
    for (s = 0; s < label->numStrata; ++s) {
      for (l = 0; l < label->stratumSizes[s]; l++) {
        ierr = PetscSectionGetDof(rootSection, label->points[s][l], &dof);CHKERRQ(ierr);
        ierr = PetscSectionSetDof(rootSection, label->points[s][l], dof+1);CHKERRQ(ierr);
      }
    }
  }
  ierr = PetscSectionSetUp(rootSection);CHKERRQ(ierr);
  /* Create a point-wise array of stratum values */
  ierr = PetscSectionGetStorageSize(rootSection, &size);CHKERRQ(ierr);
  ierr = PetscMalloc1(size, &rootStrata);CHKERRQ(ierr);
  ierr = PetscCalloc1(nroots, &rootIdx);CHKERRQ(ierr);
  if (label) {
    for (s = 0; s < label->numStrata; ++s) {
      for (l = 0; l < label->stratumSizes[s]; l++) {
        const PetscInt p = label->points[s][l];
        ierr = PetscSectionGetOffset(rootSection, p, &offset);CHKERRQ(ierr);
        rootStrata[offset+rootIdx[p]++] = label->stratumValues[s];
      }
    }
  }
  /* Build SF that maps label points to remote processes */
  ierr = PetscSectionCreate(comm, leafSection);CHKERRQ(ierr);
  ierr = PetscSFDistributeSection(sf, rootSection, &remoteOffsets, *leafSection);CHKERRQ(ierr);
  ierr = PetscSFCreateSectionSF(sf, rootSection, remoteOffsets, *leafSection, &labelSF);CHKERRQ(ierr);
  ierr = PetscFree(remoteOffsets);CHKERRQ(ierr);
  /* Send the strata for each point over the derived SF */
  ierr = PetscSectionGetStorageSize(*leafSection, &size);CHKERRQ(ierr);
  ierr = PetscMalloc1(size, leafStrata);CHKERRQ(ierr);
  ierr = PetscSFBcastBegin(labelSF, MPIU_INT, rootStrata, *leafStrata);CHKERRQ(ierr);
  ierr = PetscSFBcastEnd(labelSF, MPIU_INT, rootStrata, *leafStrata);CHKERRQ(ierr);
  /* Clean up */
  ierr = PetscFree(rootStrata);CHKERRQ(ierr);
  ierr = PetscFree(rootIdx);CHKERRQ(ierr);
  ierr = PetscSectionDestroy(&rootSection);CHKERRQ(ierr);
  ierr = PetscSFDestroy(&labelSF);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*@C
   PetscViewersCreate - Creates a container to hold a set of PetscViewers.

   Collective on MPI_Comm

   Input Parameter:
.   comm - the MPI communicator

   Output Parameter:
.  v - the collection of PetscViewers

   Level: intermediate

   Concepts: PetscViewer^array of

.seealso: PetscViewerCreate(), PetscViewersDestroy()

@*/
PetscErrorCode  PetscViewersCreate(MPI_Comm comm,PetscViewers *v)
{
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr       = PetscNew(v);CHKERRQ(ierr);
  (*v)->n    = 64;
  (*v)->comm = comm;

  ierr = PetscCalloc1(64,&(*v)->viewer);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PETSC_EXTERN PetscErrorCode PetscDrawCreate_Image(PetscDraw draw)
{
  PetscImage     img;
  int            w = draw->w, h = draw->h;
  PetscInt       size[2], nsize = 2;
  PetscBool      set;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  draw->pause   = 0;
  draw->coor_xl = 0; draw->coor_xr = 1;
  draw->coor_yl = 0; draw->coor_yr = 1;
  draw->port_xl = 0; draw->port_xr = 1;
  draw->port_yl = 0; draw->port_yr = 1;

  size[0] = w; if (size[0] < 1) size[0] = 300;
  size[1] = h; if (size[1] < 1) size[1] = size[0];
  ierr = PetscOptionsGetIntArray(((PetscObject)draw)->options,((PetscObject)draw)->prefix,"-draw_size",size,&nsize,&set);CHKERRQ(ierr);
  if (set && nsize == 1) size[1] = size[0];
  if (size[0] < 1) size[0] = 300;
  if (size[1] < 1) size[1] = size[0];
  draw->w = w = size[0]; draw->x = 0;
  draw->h = h = size[1]; draw->x = 0;

  ierr = PetscNewLog(draw,&img);CHKERRQ(ierr);
  ierr = PetscMemcpy(draw->ops,&DvOps,sizeof(DvOps));CHKERRQ(ierr);
  draw->data = (void*)img;

  img->w = w; img->h = h;
  ierr = PetscCalloc1((size_t)(img->w*img->h),&img->buffer);CHKERRQ(ierr);
  PetscImageSetClip(img,0,0,img->w,img->h);
  {
    int i,k,ncolors = 256-PETSC_DRAW_BASIC_COLORS;
    unsigned char R[256-PETSC_DRAW_BASIC_COLORS];
    unsigned char G[256-PETSC_DRAW_BASIC_COLORS];
    unsigned char B[256-PETSC_DRAW_BASIC_COLORS];
    ierr = PetscDrawUtilitySetCmap(NULL,ncolors,R,G,B);CHKERRQ(ierr);
    for (k=0; k<PETSC_DRAW_BASIC_COLORS; k++) {
      img->palette[k][0] = BasicColors[k][0];
      img->palette[k][1] = BasicColors[k][1];
      img->palette[k][2] = BasicColors[k][2];
    }
    for (i=0; i<ncolors; i++, k++) {
      img->palette[k][0] = R[i];
      img->palette[k][1] = G[i];
      img->palette[k][2] = B[i];
    }
  }

  if (!draw->savefilename) {ierr = PetscDrawSetSave(draw,draw->title);CHKERRQ(ierr);}
  PetscFunctionReturn(0);
}

/*@C
  PetscSubcommSetTypeGeneral - Set a PetscSubcomm from user's specifications

   Collective on MPI_Comm

   Input Parameter:
+  psubcomm - PetscSubcomm context
.  color   - control of subset assignment (nonnegative integer). Processes with the same color are in the same subcommunicator.
-  subrank - rank in the subcommunicator

   Level: advanced

.keywords: communicator, create

.seealso: PetscSubcommCreate(),PetscSubcommDestroy(),PetscSubcommSetNumber(),PetscSubcommSetType()
@*/
PetscErrorCode PetscSubcommSetTypeGeneral(PetscSubcomm psubcomm,PetscMPIInt color,PetscMPIInt subrank)
{
  PetscErrorCode ierr;
  MPI_Comm       subcomm=0,dupcomm=0,comm=psubcomm->parent;
  PetscMPIInt    size,icolor,duprank,*recvbuf,sendbuf[3],mysubsize,rank,*subsize;
  PetscMPIInt    i,nsubcomm=psubcomm->n;

  PetscFunctionBegin;
  if (!psubcomm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"PetscSubcomm is not created. Call PetscSubcommCreate()");
  if (nsubcomm < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"number of subcommunicators %d is incorrect. Call PetscSubcommSetNumber()",nsubcomm);

  ierr = MPI_Comm_split(comm,color,subrank,&subcomm);CHKERRQ(ierr);

  /* create dupcomm with same size as comm, but its rank, duprank, maps subcomm's contiguously into dupcomm */
  /* TODO: this can be done in an ostensibly scalale way (i.e., without allocating an array of size 'size') as is done in PetscObjectsCreateGlobalOrdering(). */
  ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
  ierr = PetscMalloc1(2*size,&recvbuf);CHKERRQ(ierr);

  ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
  ierr = MPI_Comm_size(subcomm,&mysubsize);CHKERRQ(ierr);

  sendbuf[0] = color;
  sendbuf[1] = mysubsize;
  ierr = MPI_Allgather(sendbuf,2,MPI_INT,recvbuf,2,MPI_INT,comm);CHKERRQ(ierr);

  ierr = PetscCalloc1(nsubcomm,&subsize);CHKERRQ(ierr);
  for (i=0; i<2*size; i+=2) {
    subsize[recvbuf[i]] = recvbuf[i+1];
  }
  ierr = PetscFree(recvbuf);CHKERRQ(ierr);

  duprank = 0;
  for (icolor=0; icolor<nsubcomm; icolor++) {
    if (icolor != color) { /* not color of this process */
      duprank += subsize[icolor];
    } else {
      duprank += subrank;
      break;
    }
  }
  ierr = MPI_Comm_split(comm,0,duprank,&dupcomm);CHKERRQ(ierr);

  ierr = PetscCommDuplicate(dupcomm,&psubcomm->dupparent,NULL);CHKERRQ(ierr);
  ierr = PetscCommDuplicate(subcomm,&psubcomm->child,NULL);CHKERRQ(ierr);
  ierr = MPI_Comm_free(&dupcomm);CHKERRQ(ierr);
  ierr = MPI_Comm_free(&subcomm);CHKERRQ(ierr);

  psubcomm->color   = color;
  psubcomm->subsize = subsize;
  psubcomm->type    = PETSC_SUBCOMM_GENERAL;
  PetscFunctionReturn(0);
}

/*@C
   ISColoringGetIS - Extracts index sets from the coloring context

   Collective on ISColoring

   Input Parameter:
.  iscoloring - the coloring context

   Output Parameters:
+  nn - number of index sets in the coloring context
-  is - array of index sets

   Level: advanced

.seealso: ISColoringRestoreIS(), ISColoringView()
@*/
PetscErrorCode  ISColoringGetIS(ISColoring iscoloring,PetscInt *nn,IS *isis[])
{
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscValidPointer(iscoloring,1);

  if (nn) *nn = iscoloring->n;
  if (isis) {
    if (!iscoloring->is) {
      PetscInt        *mcolors,**ii,nc = iscoloring->n,i,base, n = iscoloring->N;
      ISColoringValue *colors = iscoloring->colors;
      IS              *is;

#if defined(PETSC_USE_DEBUG)
      for (i=0; i<n; i++) {
        if (((PetscInt)colors[i]) >= nc) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coloring is our of range index %d value %d number colors %d",(int)i,(int)colors[i],(int)nc);
      }
#endif

      /* generate the lists of nodes for each color */
      ierr = PetscCalloc1(nc,&mcolors);CHKERRQ(ierr);
      for (i=0; i<n; i++) mcolors[colors[i]]++;

      ierr = PetscMalloc1(nc,&ii);CHKERRQ(ierr);
      ierr = PetscMalloc1(n,&ii[0]);CHKERRQ(ierr);
      for (i=1; i<nc; i++) ii[i] = ii[i-1] + mcolors[i-1];
      ierr = PetscMemzero(mcolors,nc*sizeof(PetscInt));CHKERRQ(ierr);

      if (iscoloring->ctype == IS_COLORING_GLOBAL) {
        ierr = MPI_Scan(&iscoloring->N,&base,1,MPIU_INT,MPI_SUM,iscoloring->comm);CHKERRQ(ierr);
        base -= iscoloring->N;
        for (i=0; i<n; i++) ii[colors[i]][mcolors[colors[i]]++] = i + base; /* global idx */
      } else if (iscoloring->ctype == IS_COLORING_GHOSTED) {
        for (i=0; i<n; i++) ii[colors[i]][mcolors[colors[i]]++] = i;   /* local idx */
      } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not provided for this ISColoringType type");

      ierr = PetscMalloc1(nc,&is);CHKERRQ(ierr);
      for (i=0; i<nc; i++) {
        ierr = ISCreateGeneral(iscoloring->comm,mcolors[i],ii[i],PETSC_COPY_VALUES,is+i);CHKERRQ(ierr);
      }

      iscoloring->is = is;
      ierr = PetscFree(ii[0]);CHKERRQ(ierr);
      ierr = PetscFree(ii);CHKERRQ(ierr);
      ierr = PetscFree(mcolors);CHKERRQ(ierr);
    }
    *isis = iscoloring->is;
  }
  PetscFunctionReturn(0);
}

static PetscErrorCode PetscDrawResizeWindow_Image(PetscDraw draw,int w,int h)
{
  PetscImage     img = (PetscImage)draw->data;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  if (w == img->w && h == img->h) PetscFunctionReturn(0);
  ierr = PetscFree(img->buffer);CHKERRQ(ierr);

  img->w = w; img->h = h;
  ierr = PetscCalloc1((size_t)(img->w*img->h),&img->buffer);CHKERRQ(ierr);
  ierr = PetscDrawSetViewport_Image(draw,draw->port_xl,draw->port_yl,draw->port_xr,draw->port_yr);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

static PetscErrorCode CreatePoints_GridReplicated(DM dm, PetscInt *Np, PetscReal **pcoords, PetscBool *pointsAllProcs, AppCtx *ctx)
{
  PetscInt       N = 3, n = 0, spaceDim, i, j, k, *ind, d;
  PetscReal      *h;
  PetscMPIInt    rank;
  PetscErrorCode ierr;

  ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRQ(ierr);
  ierr = DMGetCoordinateDim(dm, &spaceDim);CHKERRQ(ierr);
  ierr = PetscCalloc1(spaceDim,&ind);CHKERRQ(ierr);
  ierr = PetscCalloc1(spaceDim,&h);CHKERRQ(ierr);
  h[0] = 1.0/(N-1); h[1] = 1.0/(N-1); h[2] = 1.0/(N-1);
  *Np  = N * (ctx->dim > 1 ? N : 1) * (ctx->dim > 2 ? N : 1);
  ierr = PetscCalloc1(*Np * spaceDim, pcoords);CHKERRQ(ierr);
  for (k = 0; k < N; ++k) {
    ind[2] = k;
    for (j = 0; j < N; ++j) {
      ind[1] = j;
      for (i = 0; i < N; ++i, ++n) {
        ind[0] = i;

        for (d = 0; d < spaceDim; ++d) (*pcoords)[n*spaceDim+d] = ind[d]*h[d];
        ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %D (", rank, n);CHKERRQ(ierr);
        for (d = 0; d < spaceDim; ++d) {
          ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "%g", (double)(*pcoords)[n*spaceDim+d]);CHKERRQ(ierr);
          if (d < spaceDim-1) {ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, ", ");CHKERRQ(ierr);}
        }
        ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, ")\n");CHKERRQ(ierr);
      }
    }
  }
  ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL);CHKERRQ(ierr);
  *pointsAllProcs = PETSC_TRUE;
  ierr = PetscFree(ind);CHKERRQ(ierr);
  ierr = PetscFree(h);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*
 MatGetColumnIJ_SeqSELL_Color() and MatRestoreColumnIJ_SeqSELL_Color() are customized from
 MatGetColumnIJ_SeqSELL() and MatRestoreColumnIJ_SeqSELL() by adding an output
 spidx[], index of a->a, to be used in MatTransposeColoringCreate_SeqSELL() and MatFDColoringCreate_SeqSELL()
*/
PetscErrorCode MatGetColumnIJ_SeqSELL_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done)
{
  Mat_SeqSELL    *a = (Mat_SeqSELL*)A->data;
  PetscInt       i,j,*collengths,*cia,*cja,n = A->cmap->n,totalslices;
  PetscInt       row,col;
  PetscInt       *cspidx;
  PetscBool      isnonzero;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  *nn = n;
  if (!ia) PetscFunctionReturn(0);

  ierr = PetscCalloc1(n+1,&collengths);CHKERRQ(ierr);
  ierr = PetscMalloc1(n+1,&cia);CHKERRQ(ierr);
  ierr = PetscMalloc1(a->nz+1,&cja);CHKERRQ(ierr);
  ierr = PetscMalloc1(a->nz+1,&cspidx);CHKERRQ(ierr);

  totalslices = A->rmap->n/8+((A->rmap->n & 0x07)?1:0); /* floor(n/8) */
  for (i=0; i<totalslices; i++) { /* loop over slices */
    for (j=a->sliidx[i],row=0; j<a->sliidx[i+1]; j++,row=((row+1)&0x07)) {
      isnonzero = (PetscBool)((j-a->sliidx[i])/8 < a->rlen[8*i+row]);
      if (isnonzero) collengths[a->colidx[j]]++;
    }
  }

  cia[0] = oshift;
  for (i=0; i<n; i++) {
    cia[i+1] = cia[i] + collengths[i];
  }
  ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr);

  for (i=0; i<totalslices; i++) { /* loop over slices */
    for (j=a->sliidx[i],row=0; j<a->sliidx[i+1]; j++,row=((row+1)&0x07)) {
      isnonzero = (PetscBool)((j-a->sliidx[i])/8 < a->rlen[8*i+row]);
      if (isnonzero) {
        col = a->colidx[j];
        cspidx[cia[col]+collengths[col]-oshift] = j; /* index of a->colidx */
        cja[cia[col]+collengths[col]-oshift] = 8*i+row +oshift; /* row index */
        collengths[col]++;
      }
    }
  }

  ierr   = PetscFree(collengths);CHKERRQ(ierr);
  *ia    = cia; *ja = cja;
  *spidx = cspidx;
  PetscFunctionReturn(0);
}

PetscErrorCode MatPartitioningSetFromOptions_Hierarchical(PetscOptions *PetscOptionsObject,MatPartitioning part)
{
  MatPartitioning_Hierarchical *hpart = (MatPartitioning_Hierarchical*)part->data;
  PetscErrorCode ierr;
  char           value[1024];
  PetscBool      flag = PETSC_FALSE;

  PetscFunctionBegin;
  ierr = PetscOptionsHead(PetscOptionsObject,"Set hierarchical partitioning options");CHKERRQ(ierr);
  ierr = PetscOptionsString("-mat_partitioning_hierarchical_coarseparttype","coarse part type",PETSC_NULL,PETSC_NULL,value,1024,&flag);CHKERRQ(ierr);
  if(flag){
   ierr = PetscCalloc1(1024,&hpart->coarseparttype);CHKERRQ(ierr);
   ierr = PetscStrcpy(hpart->coarseparttype,value);CHKERRQ(ierr);
  }
  ierr = PetscOptionsString("-mat_partitioning_hierarchical_fineparttype","fine part type",PETSC_NULL,PETSC_NULL,value,1024,&flag);CHKERRQ(ierr);
  if(flag){
    ierr = PetscCalloc1(1024,&hpart->fineparttype);CHKERRQ(ierr);
    ierr = PetscStrcpy(hpart->fineparttype,value);CHKERRQ(ierr);
  }
  ierr = PetscOptionsInt("-mat_partitioning_hierarchical_Ncoarseparts","number of coarse parts",PETSC_NULL,0,&hpart->Ncoarseparts,&flag);CHKERRQ(ierr);
  ierr = PetscOptionsInt("-mat_partitioning_hierarchical_Nfineparts","number of fine parts",PETSC_NULL,1,&hpart->Nfineparts,&flag);CHKERRQ(ierr);
  ierr = PetscOptionsTail();CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

static PetscErrorCode PCSetUp_CP(PC pc)
{
  PC_CP          *cp = (PC_CP*)pc->data;
  PetscInt       i,j,*colcnt;
  PetscErrorCode ierr;
  PetscBool      flg;
  Mat_SeqAIJ     *aij = (Mat_SeqAIJ*)pc->pmat->data;

  PetscFunctionBegin;
  ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATSEQAIJ,&flg);CHKERRQ(ierr);
  if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Currently only handles SeqAIJ matrices");

  ierr = MatGetLocalSize(pc->pmat,&cp->m,&cp->n);CHKERRQ(ierr);
  if (cp->m != cp->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently only for square matrices");

  if (!cp->work) {ierr = MatCreateVecs(pc->pmat,&cp->work,NULL);CHKERRQ(ierr);}
  if (!cp->d) {ierr = PetscMalloc1(cp->n,&cp->d);CHKERRQ(ierr);}
  if (cp->a && pc->flag != SAME_NONZERO_PATTERN) {
    ierr  = PetscFree3(cp->a,cp->i,cp->j);CHKERRQ(ierr);
    cp->a = 0;
  }

  /* convert to column format */
  if (!cp->a) {
    ierr = PetscMalloc3(aij->nz,&cp->a,cp->n+1,&cp->i,aij->nz,&cp->j);CHKERRQ(ierr);
  }
  ierr = PetscCalloc1(cp->n,&colcnt);CHKERRQ(ierr);

  for (i=0; i<aij->nz; i++) colcnt[aij->j[i]]++;
  cp->i[0] = 0;
  for (i=0; i<cp->n; i++) cp->i[i+1] = cp->i[i] + colcnt[i];
  ierr = PetscMemzero(colcnt,cp->n*sizeof(PetscInt));CHKERRQ(ierr);
  for (i=0; i<cp->m; i++) {  /* over rows */
    for (j=aij->i[i]; j<aij->i[i+1]; j++) {  /* over columns in row */
      cp->j[cp->i[aij->j[j]]+colcnt[aij->j[j]]]   = i;
      cp->a[cp->i[aij->j[j]]+colcnt[aij->j[j]]++] = aij->a[j];
    }
  }
  ierr = PetscFree(colcnt);CHKERRQ(ierr);

  /* compute sum of squares of each column d[] */
  for (i=0; i<cp->n; i++) {  /* over columns */
    cp->d[i] = 0.;
    for (j=cp->i[i]; j<cp->i[i+1]; j++) cp->d[i] += cp->a[j]*cp->a[j]; /* over rows in column */
    cp->d[i] = 1.0/cp->d[i];
  }
  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 = PetscCalloc1((newmax-PETSC_SMALLEST_CLASSID),&newbase);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);
}

/*
   PetscTableCreate  Creates a PETSc look up table

   Input Parameters:
+     n - expected number of keys
-     maxkey- largest possible key

   Notes: keys are between 1 and maxkey inclusive

*/
PetscErrorCode  PetscTableCreate(const PetscInt n,PetscInt maxkey,PetscTable *rta)
{
  PetscTable     ta;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  if (n < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"n < 0");
  ierr          = PetscNew(&ta);CHKERRQ(ierr);
  ta->tablesize = (3*n)/2 + 17;
  if (ta->tablesize < n) ta->tablesize = PETSC_MAX_INT/4; /* overflow */
  ierr       = PetscCalloc1(ta->tablesize,&ta->keytable);CHKERRQ(ierr);
  ierr       = PetscMalloc1(ta->tablesize,&ta->table);CHKERRQ(ierr);
  ta->head   = 0;
  ta->count  = 0;
  ta->maxkey = maxkey;
  *rta       = ta;
  PetscFunctionReturn(0);
}

/*@C
  DMPlexAddLabel - Add the label to this mesh

  Not Collective

  Input Parameters:
+ dm   - The DMPlex object
- label - The DMLabel

  Level: developer

.keywords: mesh
.seealso: DMPlexCreateLabel(), DMPlexHasLabel(), DMPlexGetLabelValue(), DMPlexSetLabelValue(), DMPlexGetStratumIS()
@*/
PetscErrorCode DMPlexAddLabel(DM dm, DMLabel label)
{
  DM_Plex       *mesh = (DM_Plex*) dm->data;
  PlexLabel      tmpLabel;
  PetscBool      hasLabel;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
  ierr = DMPlexHasLabel(dm, label->name, &hasLabel);CHKERRQ(ierr);
  if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", label->name);
  ierr = PetscCalloc1(1, &tmpLabel);CHKERRQ(ierr);
  tmpLabel->label  = label;
  tmpLabel->output = PETSC_TRUE;
  tmpLabel->next   = mesh->labels;
  mesh->labels     = tmpLabel;
  PetscFunctionReturn(0);
}