/*@
   DMRestoreLocalVector - Returns a Seq PETSc vector that
     obtained from DMGetLocalVector(). Do not use with vector obtained via
     DMCreateLocalVector().

   Not Collective

   Input Parameter:
+  dm - the distributed array
-  g - the local vector

   Level: beginner

.keywords: distributed array, create, local, vector

.seealso: DMCreateGlobalVector(), VecDuplicate(), VecDuplicateVecs(),
          DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMGlobalToLocalBegin(),
          DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMCreateLocalVector(), DMGetLocalVector()
@*/
PetscErrorCode  DMRestoreLocalVector(DM dm,Vec* g)
{
  PetscErrorCode ierr;
  PetscInt       i,j;

  PetscFunctionBegin; 
  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
  PetscValidPointer(g,2);
  for (j=0; j<DM_MAX_WORK_VECTORS; j++) {
    if (*g == dm->localout[j]) {
      dm->localout[j] = PETSC_NULL;
      for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
        if (!dm->localin[i]) {
          dm->localin[i] = *g;
          goto alldone;
        }
      }
    }
  }
  ierr = VecDestroy(g);CHKERRQ(ierr);
  alldone:
  PetscFunctionReturn(0);
}

/*@C
   PetscDrawLGSetLegend - sets the names of each curve plotted

   Logically Collective over PetscDrawLG

   Input Parameter:
+  lg - the line graph context.
-  names - the names for each curve

   Level: intermediate

   Concepts: line graph^setting number of lines

@*/
PetscErrorCode  PetscDrawLGSetLegend(PetscDrawLG lg,const char *const *names)
{
  PetscErrorCode ierr;
  PetscInt       i;

  PetscFunctionBegin;
  if (lg && ((PetscObject)lg)->classid == PETSC_DRAW_CLASSID) PetscFunctionReturn(0);
  PetscValidHeaderSpecific(lg,PETSC_DRAWLG_CLASSID,1);

  if (lg->legend) {
    for (i=0; i<lg->dim; i++) {
      ierr = PetscFree(lg->legend[i]);CHKERRQ(ierr);
    }
    ierr = PetscFree(lg->legend);CHKERRQ(ierr);
  }
  if (names) {
    ierr = PetscMalloc(lg->dim*sizeof(char**),&lg->legend);CHKERRQ(ierr);
    for (i=0; i<lg->dim; i++) {
      ierr = PetscStrallocpy(names[i],&lg->legend[i]);CHKERRQ(ierr);
    }
  }
  PetscFunctionReturn(0);
}

/*@C
   PetscSFSetType - set the PetscSF communication implementation

   Collective on PetscSF

   Input Parameters:
+  sf - the PetscSF context
-  type - a known method

   Options Database Key:
.  -sf_type <type> - Sets the method; use -help for a list
   of available methods (for instance, window, pt2pt, neighbor)

   Notes:
   See "include/petscsf.h" for available methods (for instance)
+    PETSCSFWINDOW - MPI-2/3 one-sided
-    PETSCSFBASIC - basic implementation using MPI-1 two-sided

  Level: intermediate

.keywords: PetscSF, set, type

.seealso: PetscSFType, PetscSFCreate()
@*/
PetscErrorCode PetscSFSetType(PetscSF sf,PetscSFType type)
{
  PetscErrorCode ierr,(*r)(PetscSF);
  PetscBool      match;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(sf,PETSCSF_CLASSID,1);
  PetscValidCharPointer(type,2);

  ierr = PetscObjectTypeCompare((PetscObject)sf,type,&match);CHKERRQ(ierr);
  if (match) PetscFunctionReturn(0);

  ierr = PetscFunctionListFind(PetscSFList,type,&r);CHKERRQ(ierr);
  if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested PetscSF type %s",type);
  /* Destroy the previous private PetscSF context */
  if (sf->ops->Destroy) {
    ierr = (*(sf)->ops->Destroy)(sf);CHKERRQ(ierr);
  }
  ierr = PetscMemzero(sf->ops,sizeof(*sf->ops));CHKERRQ(ierr);
  ierr = PetscObjectChangeTypeName((PetscObject)sf,type);CHKERRQ(ierr);
  ierr = (*r)(sf);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*@
   MatFDColoringSetFromOptions - Sets coloring finite difference parameters from
   the options database.

   Collective on MatFDColoring

   The Jacobian, F'(u), is estimated with the differencing approximation
.vb
       F'(u)_{:,i} = [F(u+h*dx_{i}) - F(u)]/h where
       h = error_rel*u[i]                 if  abs(u[i]) > umin
         = +/- error_rel*umin             otherwise, with +/- determined by the sign of u[i]
       dx_{i} = (0, ... 1, .... 0)
.ve

   Input Parameter:
.  coloring - the coloring context

   Options Database Keys:
+  -mat_fd_coloring_err <err> - Sets <err> (square root
           of relative error in the function)
.  -mat_fd_coloring_umin <umin> - Sets umin, the minimum allowable u-value magnitude
.  -mat_fd_type - "wp" or "ds" (see MATMFFD_WP or MATMFFD_DS)
.  -mat_fd_coloring_view - Activates basic viewing
.  -mat_fd_coloring_view ::ascii_info - Activates viewing info
-  -mat_fd_coloring_view draw - Activates drawing

    Level: intermediate

.keywords: Mat, finite differences, parameters

.seealso: MatFDColoringCreate(), MatFDColoringView(), MatFDColoringSetParameters()

@*/
PetscErrorCode  MatFDColoringSetFromOptions(MatFDColoring matfd)
{
  PetscErrorCode ierr;
  PetscBool      flg;
  char           value[3];

  PetscFunctionBegin;
  PetscValidHeaderSpecific(matfd,MAT_FDCOLORING_CLASSID,1);

  ierr = PetscObjectOptionsBegin((PetscObject)matfd);CHKERRQ(ierr);
    ierr = PetscOptionsReal("-mat_fd_coloring_err","Square root of relative error in function","MatFDColoringSetParameters",matfd->error_rel,&matfd->error_rel,0);CHKERRQ(ierr);
    ierr = PetscOptionsReal("-mat_fd_coloring_umin","Minimum allowable u magnitude","MatFDColoringSetParameters",matfd->umin,&matfd->umin,0);CHKERRQ(ierr);
    ierr = PetscOptionsString("-mat_fd_type","Algorithm to compute h, wp or ds","MatFDColoringCreate",matfd->htype,value,3,&flg);CHKERRQ(ierr);
    if (flg) {
      if (value[0] == 'w' && value[1] == 'p') matfd->htype = "wp";
      else if (value[0] == 'd' && value[1] == 's') matfd->htype = "ds";
      else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Unknown finite differencing type %s",value);
    }
    /* process any options handlers added with PetscObjectAddOptionsHandler() */
    ierr = PetscObjectProcessOptionsHandlers((PetscObject)matfd);CHKERRQ(ierr);
  PetscOptionsEnd();CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode DMCreateLocalVector_Shell(DM dm,Vec *gvec)
{
  PetscErrorCode ierr;
  DM_Shell       *shell = (DM_Shell*)dm->data;
  Vec            X;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
  PetscValidPointer(gvec,2);
  *gvec = 0;
  X     = shell->Xlocal;
  if (!X) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMShellSetLocalVector() or DMShellSetCreateLocalVector()");
  if (((PetscObject)X)->refct < 2) { /* We have an exclusive reference so we can give it out */
    ierr  = PetscObjectReference((PetscObject)X);CHKERRQ(ierr);
    ierr  = VecZeroEntries(X);CHKERRQ(ierr);
    *gvec = X;
  } else {                      /* Need to create a copy, could use MAT_SHARE_NONZERO_PATTERN in most cases */
    ierr = VecDuplicate(X,gvec);CHKERRQ(ierr);
    ierr = VecZeroEntries(*gvec);CHKERRQ(ierr);
  }
  ierr = VecSetDM(*gvec,dm);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode MatCoarsenView_MIS(MatCoarsen coarse,PetscViewer viewer)
{
    /* MatCoarsen_MIS *MIS = (MatCoarsen_MIS*)coarse->; */
    PetscErrorCode ierr;
    int            rank;
    PetscBool      iascii;

    PetscFunctionBegin;
    PetscValidHeaderSpecific(coarse,MAT_COARSEN_CLASSID,1);
    ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)coarse),&rank);
    CHKERRQ(ierr);
    ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
    CHKERRQ(ierr);
    if (iascii) {
        ierr = PetscViewerASCIISynchronizedPrintf(viewer,"  [%d] MIS aggregator\n",rank);
        CHKERRQ(ierr);
        ierr = PetscViewerFlush(viewer);
        CHKERRQ(ierr);
        ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
        CHKERRQ(ierr);
    }
    PetscFunctionReturn(0);
}

dErr dUnitsDestroy(dUnits *unp)
{
  dUnits un = *unp;
  dErr err;

  dFunctionBegin;
  if (!un) dFunctionReturn(0);
  PetscValidHeaderSpecific(un,dUNITS_CLASSID,1);
  if (--((PetscObject)un)->refct > 0) dFunctionReturn(0);
  for (dInt i=0; i<un->nalloc; i++) {
    dUnit u = un->list[i];
    if (u) {
      err = dFree(u->quantity);dCHK(err);
      err = dFree(u->longname);dCHK(err);
      err = dFree(u->shortname);dCHK(err);
      err = dFree(u->siname);dCHK(err);
    }
    err = dFree(un->list[i]);dCHK(err);
  }
  err = dFree(un->list);dCHK(err);
  err = PetscHeaderDestroy(unp);dCHK(err);
  dFunctionReturn(0);
}

/*@
  TaoLineSearchSetFromOptions - Sets various TaoLineSearch parameters from user
  options.

  Collective on TaoLineSearch

  Input Paremeter:
. ls - the TaoLineSearch context

  Options Database Keys:
+ -tao_ls_type <type> - The algorithm that TAO uses (more-thuente, gpcg, unit)
. -tao_ls_ftol <tol> - tolerance for sufficient decrease
. -tao_ls_gtol <tol> - tolerance for curvature condition
. -tao_ls_rtol <tol> - relative tolerance for acceptable step
. -tao_ls_stepmin <step> - minimum steplength allowed
. -tao_ls_stepmax <step> - maximum steplength allowed
. -tao_ls_max_funcs <n> - maximum number of function evaluations allowed
- -tao_ls_view - display line-search results to standard output

  Level: beginner
@*/
PetscErrorCode TaoLineSearchSetFromOptions(TaoLineSearch ls)
{
  PetscErrorCode ierr;
  const char     *default_type=TAOLINESEARCHMT;
  char           type[256];
  PetscBool      flg;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ls,TAOLINESEARCH_CLASSID,1);
  ierr = PetscObjectOptionsBegin((PetscObject)ls);CHKERRQ(ierr);
  if (!TaoLineSearchInitialized) {
    ierr = TaoLineSearchInitializePackage();CHKERRQ(ierr);
  }
  if (((PetscObject)ls)->type_name) {
    default_type = ((PetscObject)ls)->type_name;
  }
  /* Check for type from options */
  ierr = PetscOptionsFList("-tao_ls_type","Tao Line Search type","TaoLineSearchSetType",TaoLineSearchList,default_type,type,256,&flg);CHKERRQ(ierr);
  if (flg) {
    ierr = TaoLineSearchSetType(ls,type);CHKERRQ(ierr);
  } else if (!((PetscObject)ls)->type_name) {
    ierr = TaoLineSearchSetType(ls,default_type);
  }

  ierr = PetscOptionsInt("-tao_ls_max_funcs","max function evals in line search","",ls->max_funcs,&ls->max_funcs,0);CHKERRQ(ierr);
  ierr = PetscOptionsReal("-tao_ls_ftol","tol for sufficient decrease","",ls->ftol,&ls->ftol,0);CHKERRQ(ierr);
  ierr = PetscOptionsReal("-tao_ls_gtol","tol for curvature condition","",ls->gtol,&ls->gtol,0);CHKERRQ(ierr);
  ierr = PetscOptionsReal("-tao_ls_rtol","relative tol for acceptable step","",ls->rtol,&ls->rtol,0);CHKERRQ(ierr);
  ierr = PetscOptionsReal("-tao_ls_stepmin","lower bound for step","",ls->stepmin,&ls->stepmin,0);CHKERRQ(ierr);
  ierr = PetscOptionsReal("-tao_ls_stepmax","upper bound for step","",ls->stepmax,&ls->stepmax,0);CHKERRQ(ierr);
  ierr = PetscOptionsBool("-tao_ls_view","view TaoLineSearch info after each line search has completed","TaoLineSearchView",PETSC_FALSE,&ls->viewls,NULL);CHKERRQ(ierr);
  if (ls->ops->setfromoptions) {
    ierr = (*ls->ops->setfromoptions)(ls);CHKERRQ(ierr);
  }
  ierr = PetscOptionsEnd();CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*@
   EPSReset - Resets the EPS context to the initial state and removes any
   allocated objects.

   Collective on EPS

   Input Parameter:
.  eps - eigensolver context obtained from EPSCreate()

   Level: advanced

.seealso: EPSDestroy()
@*/
PetscErrorCode EPSReset(EPS eps)
{
    PetscErrorCode ierr;
    PetscInt       ncols;

    PetscFunctionBegin;
    PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
    if (eps->ops->reset) {
        ierr = (eps->ops->reset)(eps);
        CHKERRQ(ierr);
    }
    if (eps->st) {
        ierr = STReset(eps->st);
        CHKERRQ(ierr);
    }
    if (eps->ds) {
        ierr = DSReset(eps->ds);
        CHKERRQ(ierr);
    }
    ierr = VecDestroy(&eps->D);
    CHKERRQ(ierr);
    ierr = BVGetSizes(eps->V,NULL,NULL,&ncols);
    CHKERRQ(ierr);
    if (ncols) {
        ierr = PetscFree4(eps->eigr,eps->eigi,eps->errest,eps->perm);
        CHKERRQ(ierr);
        ierr = PetscFree2(eps->rr,eps->ri);
        CHKERRQ(ierr);
    }
    ierr = BVDestroy(&eps->V);
    CHKERRQ(ierr);
    ierr = VecDestroyVecs(eps->nwork,&eps->work);
    CHKERRQ(ierr);
    eps->nwork = 0;
    eps->state = EPS_STATE_INITIAL;
    PetscFunctionReturn(0);
}

/*@
   PCSetFromOptions - Sets PC options from the options database.
   This routine must be called before PCSetUp() if the user is to be
   allowed to set the preconditioner method.

   Collective on PC

   Input Parameter:
.  pc - the preconditioner context

   Options Database:
.   -pc_use_amat true,false see PCSetUseAmat()

   Level: developer

.keywords: PC, set, from, options, database

.seealso: PCSetUseAmat()

@*/
PetscErrorCode  PCSetFromOptions(PC pc)
{
  PetscErrorCode ierr;
  char           type[256];
  const char     *def;
  PetscBool      flg;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc,PC_CLASSID,1);

  if (!PCRegisterAllCalled) {ierr = PCRegisterAll();CHKERRQ(ierr);}
  ierr = PetscObjectOptionsBegin((PetscObject)pc);CHKERRQ(ierr);
  if (!((PetscObject)pc)->type_name) {
    ierr = PCGetDefaultType_Private(pc,&def);CHKERRQ(ierr);
  } else {
    def = ((PetscObject)pc)->type_name;
  }

  ierr = PetscOptionsList("-pc_type","Preconditioner","PCSetType",PCList,def,type,256,&flg);CHKERRQ(ierr);
  if (flg) {
    ierr = PCSetType(pc,type);CHKERRQ(ierr);
  } else if (!((PetscObject)pc)->type_name) {
    ierr = PCSetType(pc,def);CHKERRQ(ierr);
  }

  ierr = PetscOptionsBool("-pc_use_amat","use Amat (instead of Pmat) to define preconditioner in nested inner solves","PCSetUseAmat",pc->useAmat,&pc->useAmat,NULL);CHKERRQ(ierr);

  if (pc->ops->setfromoptions) {
    ierr = (*pc->ops->setfromoptions)(pc);CHKERRQ(ierr);
  }

  /* process any options handlers added with PetscObjectAddOptionsHandler() */
  ierr = PetscObjectProcessOptionsHandlers((PetscObject)pc);CHKERRQ(ierr);
  ierr = PetscOptionsEnd();CHKERRQ(ierr);
  pc->setfromoptionscalled++;
  PetscFunctionReturn(0);
}

/*@
   MatNullSpaceCreate - Creates a data structure used to project vectors
   out of null spaces.

   Collective on MPI_Comm

   Input Parameters:
+  comm - the MPI communicator associated with the object
.  has_cnst - PETSC_TRUE if the null space contains the constant vector; otherwise PETSC_FALSE
.  n - number of vectors (excluding constant vector) in null space
-  vecs - the vectors that span the null space (excluding the constant vector);
          these vectors must be orthonormal. These vectors are NOT copied, so do not change them
          after this call. You should free the array that you pass in and destroy the vectors (this will reduce the reference count
          for them by one).

   Output Parameter:
.  SP - the null space context

   Level: advanced

   Notes: See MatNullSpaceSetFunction() as an alternative way of providing the null space information instead of setting vecs.

      If has_cnst is PETSC_TRUE you do not need to pass a constant vector in as a fourth argument to this routine, nor do you
       need to pass in a function that eliminates the constant function into MatNullSpaceSetFunction().

  Users manual sections:
.   sec_singular

.keywords: PC, null space, create

.seealso: MatNullSpaceDestroy(), MatNullSpaceRemove(), KSPSetNullSpace(), MatNullSpace, MatNullSpaceSetFunction()
@*/
PetscErrorCode  MatNullSpaceCreate(MPI_Comm comm,PetscBool has_cnst,PetscInt n,const Vec vecs[],MatNullSpace *SP)
{
  MatNullSpace   sp;
  PetscErrorCode ierr;
  PetscInt       i;

  PetscFunctionBegin;
  if (n < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",n);
  if (n) PetscValidPointer(vecs,4);
  for (i=0; i<n; i++) PetscValidHeaderSpecific(vecs[i],VEC_CLASSID,4);
  PetscValidPointer(SP,5);

  *SP = NULL;
  ierr = MatInitializePackage();CHKERRQ(ierr);

  ierr = PetscHeaderCreate(sp,_p_MatNullSpace,int,MAT_NULLSPACE_CLASSID,"MatNullSpace","Null space","Mat",comm,MatNullSpaceDestroy,MatNullSpaceView);CHKERRQ(ierr);

  sp->has_cnst = has_cnst;
  sp->n        = n;
  sp->vecs     = 0;
  sp->alpha    = 0;
  sp->remove   = 0;
  sp->rmctx    = 0;

  if (n) {
    ierr = PetscMalloc1(n,&sp->vecs);CHKERRQ(ierr);
    ierr = PetscMalloc1(n,&sp->alpha);CHKERRQ(ierr);
    ierr = PetscLogObjectMemory((PetscObject)sp,n*(sizeof(Vec)+sizeof(PetscScalar)));CHKERRQ(ierr);
    for (i=0; i<n; i++) {
      ierr        = PetscObjectReference((PetscObject)vecs[i]);CHKERRQ(ierr);
      sp->vecs[i] = vecs[i];
    }
  }

  *SP = sp;
  PetscFunctionReturn(0);
}

PetscErrorCode  KSPFischerGuessUpdate_Method2(KSPFischerGuess_Method2 *itg,Vec x)
{
  PetscScalar    norm;
  PetscErrorCode ierr;
  int            curl = itg->curl,i;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(x,VEC_CLASSID,2);
  PetscValidPointer(itg,3);
  if (curl == itg->maxl) {
    ierr      = KSP_MatMult(itg->ksp,itg->mat,x,itg->Ax);CHKERRQ(ierr); /* norm = sqrt(x'Ax) */
    ierr      = VecDot(x,itg->Ax,&norm);CHKERRQ(ierr);
    ierr      = VecCopy(x,itg->xtilde[0]);CHKERRQ(ierr);
    ierr      = VecScale(itg->xtilde[0],1.0/PetscSqrtScalar(norm));CHKERRQ(ierr);
    itg->curl = 1;
  } else {
    if (!curl) {
      ierr = VecCopy(x,itg->xtilde[curl]);CHKERRQ(ierr);
    } else {
      ierr = VecWAXPY(itg->xtilde[curl],-1.0,itg->guess,x);CHKERRQ(ierr);
    }
    ierr = KSP_MatMult(itg->ksp,itg->mat,itg->xtilde[curl],itg->Ax);CHKERRQ(ierr);
    ierr = VecMDot(itg->Ax,curl,itg->xtilde,itg->alpha);CHKERRQ(ierr);
    for (i=0; i<curl; i++) itg->alpha[i] = -itg->alpha[i];
    ierr = VecMAXPY(itg->xtilde[curl],curl,itg->alpha,itg->xtilde);CHKERRQ(ierr);

    ierr = KSP_MatMult(itg->ksp,itg->mat,itg->xtilde[curl],itg->Ax);CHKERRQ(ierr); /* norm = sqrt(xtilde[curl]'Axtilde[curl]) */
    ierr = VecDot(itg->xtilde[curl],itg->Ax,&norm);CHKERRQ(ierr);
    if (PetscAbsScalar(norm) != 0.0) {
      ierr = VecScale(itg->xtilde[curl],1.0/PetscSqrtScalar(norm));CHKERRQ(ierr);
      itg->curl++;
    } else {
      ierr = PetscInfo(itg->ksp,"Not increasing dimension of Fischer space because new direction is identical to previous\n");CHKERRQ(ierr);
    }
  }
  PetscFunctionReturn(0);
}

/*@C
   PetscSFCreateEmbeddedSF - removes edges from all but the selected roots, does not remap indices

   Collective

   Input Arguments:
+  sf - original star forest
.  nroots - number of roots to select on this process
-  selected - selected roots on this process

   Output Arguments:
.  newsf - new star forest

   Level: advanced

   Note:
   To use the new PetscSF, it may be necessary to know the indices of the leaves that are still participating. This can
   be done by calling PetscSFGetGraph().

.seealso: PetscSFSetGraph(), PetscSFGetGraph()
@*/
PetscErrorCode PetscSFCreateEmbeddedSF(PetscSF sf,PetscInt nroots,const PetscInt *selected,PetscSF *newsf)
{
  PetscInt      *rootdata, *leafdata, *ilocal;
  PetscSFNode   *iremote;
  PetscInt       leafsize = 0, nleaves = 0, n, i;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(sf,PETSCSF_CLASSID,1);
  if (nroots) PetscValidPointer(selected,3);
  PetscValidPointer(newsf,4);
  if (sf->mine) for (i = 0; i < sf->nleaves; ++i) {leafsize = PetscMax(leafsize, sf->mine[i]+1);}
  else leafsize = sf->nleaves;
  ierr = PetscCalloc2(sf->nroots,&rootdata,leafsize,&leafdata);CHKERRQ(ierr);
  for (i=0; i<nroots; ++i) rootdata[selected[i]] = 1;
  ierr = PetscSFBcastBegin(sf,MPIU_INT,rootdata,leafdata);CHKERRQ(ierr);
  ierr = PetscSFBcastEnd(sf,MPIU_INT,rootdata,leafdata);CHKERRQ(ierr);

  for (i = 0; i < leafsize; ++i) nleaves += leafdata[i];
  ierr = PetscMalloc1(nleaves,&ilocal);CHKERRQ(ierr);
  ierr = PetscMalloc1(nleaves,&iremote);CHKERRQ(ierr);
  for (i = 0, n = 0; i < sf->nleaves; ++i) {
    const PetscInt lidx = sf->mine ? sf->mine[i] : i;

    if (leafdata[lidx]) {
      ilocal[n]        = lidx;
      iremote[n].rank  = sf->remote[i].rank;
      iremote[n].index = sf->remote[i].index;
      ++n;
    }
  }
  if (n != nleaves) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "There is a size mismatch in the SF embedding, %d != %d", n, nleaves);
  ierr = PetscSFDuplicate(sf,PETSCSF_DUPLICATE_RANKS,newsf);CHKERRQ(ierr);
  ierr = PetscSFSetGraph(*newsf,sf->nroots,nleaves,ilocal,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);CHKERRQ(ierr);
  ierr = PetscFree2(rootdata,leafdata);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*@C
  PetscSFCreateEmbeddedLeafSF - removes edges from all but the selected leaves, does not remap indices

  Collective

  Input Arguments:
+ sf - original star forest
. nleaves - number of leaves to select on this process
- selected - selected leaves on this process

  Output Arguments:
.  newsf - new star forest

  Level: advanced

.seealso: PetscSFCreateEmbeddedSF(), PetscSFSetGraph(), PetscSFGetGraph()
@*/
PetscErrorCode PetscSFCreateEmbeddedLeafSF(PetscSF sf, PetscInt nleaves, const PetscInt *selected, PetscSF *newsf)
{
  PetscSFNode   *iremote;
  PetscInt      *ilocal;
  PetscInt       i;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(sf, PETSCSF_CLASSID, 1);
  if (nleaves) PetscValidPointer(selected, 3);
  PetscValidPointer(newsf, 4);
  ierr = PetscMalloc1(nleaves, &ilocal);CHKERRQ(ierr);
  ierr = PetscMalloc1(nleaves, &iremote);CHKERRQ(ierr);
  for (i = 0; i < nleaves; ++i) {
    const PetscInt l = selected[i];

    ilocal[i]        = sf->mine ? sf->mine[l] : l;
    iremote[i].rank  = sf->remote[l].rank;
    iremote[i].index = sf->remote[l].index;
  }
  ierr = PetscSFDuplicate(sf, PETSCSF_DUPLICATE_RANKS, newsf);CHKERRQ(ierr);
  ierr = PetscSFSetGraph(*newsf, sf->nroots, nleaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*@C
  MatSchurComplementGetSubMatrices - Get the individual submatrices in the Schur complement

  Collective on Mat

  Input Parameter:
. S                - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01

  Output Paramters:
+ A00,A01,A10,A11  - the four parts of the original matrix A = [A00 A01; A10 A11] (A11 is optional)
- Ap00             - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}.

  Note: A11 is optional, and thus can be NULL.  The submatrices are not increfed before they are returned and should not be modified or destroyed.

  Level: intermediate

.seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatCreateSchurComplement(), MatSchurComplementUpdateSubMatrices()
@*/
PetscErrorCode  MatSchurComplementGetSubMatrices(Mat S,Mat *A00,Mat *Ap00,Mat *A01,Mat *A10,Mat *A11)
{
  Mat_SchurComplement *Na = (Mat_SchurComplement*) S->data;
  PetscErrorCode      ierr;
  PetscBool           flg;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(S,MAT_CLASSID,1);
  ierr = PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&flg);CHKERRQ(ierr);
  if (flg) {
    if (A00) *A00 = Na->A;
    if (Ap00) *Ap00 = Na->Ap;
    if (A01) *A01 = Na->B;
    if (A10) *A10 = Na->C;
    if (A11) *A11 = Na->D;
  } else {
    if (A00) *A00 = 0;
    if (Ap00) *Ap00 = 0;
    if (A01) *A01 = 0;
    if (A10) *A10 = 0;
    if (A11) *A11 = 0;
  }
  PetscFunctionReturn(0);
}

/*@C
   SNESDefaultConverged - Convergence test of the solvers for
   systems of nonlinear equations (default).

   Collective on SNES

   Input Parameters:
+  snes - the SNES context
.  it - the iteration (0 indicates before any Newton steps)
.  xnorm - 2-norm of current iterate
.  snorm - 2-norm of current step
.  fnorm - 2-norm of function at current iterate
-  dummy - unused context

   Output Parameter:
.   reason  - one of
$  SNES_CONVERGED_FNORM_ABS       - (fnorm < abstol),
$  SNES_CONVERGED_SNORM_RELATIVE  - (snorm < stol*xnorm),
$  SNES_CONVERGED_FNORM_RELATIVE  - (fnorm < rtol*fnorm0),
$  SNES_DIVERGED_FUNCTION_COUNT   - (nfct > maxf),
$  SNES_DIVERGED_FNORM_NAN        - (fnorm == NaN),
$  SNES_CONVERGED_ITERATING       - (otherwise),

   where
+    maxf - maximum number of function evaluations,
            set with SNESSetTolerances()
.    nfct - number of function evaluations,
.    abstol - absolute function norm tolerance,
            set with SNESSetTolerances()
-    rtol - relative function norm tolerance, set with SNESSetTolerances()

   Level: intermediate

.keywords: SNES, nonlinear, default, converged, convergence

.seealso: SNESSetConvergenceTest()
@*/
PetscErrorCode  SNESDefaultConverged(SNES snes,PetscInt it,PetscReal xnorm,PetscReal snorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy)
{
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
  PetscValidPointer(reason,6);

  *reason = SNES_CONVERGED_ITERATING;

  if (!it) {
    /* set parameter for default relative tolerance convergence test */
    snes->ttol = fnorm*snes->rtol;
  }
  if (PetscIsInfOrNanReal(fnorm)) {
    ierr = PetscInfo(snes,"Failed to converged, function norm is NaN\n");CHKERRQ(ierr);
    *reason = SNES_DIVERGED_FNORM_NAN;
  } else if (fnorm < snes->abstol) {
    ierr = PetscInfo2(snes,"Converged due to function norm %14.12e < %14.12e\n",(double)fnorm,(double)snes->abstol);CHKERRQ(ierr);
    *reason = SNES_CONVERGED_FNORM_ABS;
  } else if (snes->nfuncs >= snes->max_funcs) {
    ierr = PetscInfo2(snes,"Exceeded maximum number of function evaluations: %D > %D\n",snes->nfuncs,snes->max_funcs);CHKERRQ(ierr);
    *reason = SNES_DIVERGED_FUNCTION_COUNT;
  }

  if (it && !*reason) {
    if (fnorm <= snes->ttol) {
      ierr = PetscInfo2(snes,"Converged due to function norm %14.12e < %14.12e (relative tolerance)\n",(double)fnorm,(double)snes->ttol);CHKERRQ(ierr);
      *reason = SNES_CONVERGED_FNORM_RELATIVE;
    } else if (snorm < snes->stol*xnorm) {
      ierr = PetscInfo3(snes,"Converged due to small update length: %14.12e < %14.12e * %14.12e\n",(double)snorm,(double)snes->stol,(double)xnorm);CHKERRQ(ierr);
      *reason = SNES_CONVERGED_SNORM_RELATIVE;
    }
  }
  PetscFunctionReturn(0);
}

/*@C
   ISColoringView - Views a coloring context.

   Collective on ISColoring

   Input Parameters:
+  iscoloring - the coloring context
-  viewer - the viewer

   Level: advanced

.seealso: ISColoringDestroy(), ISColoringGetIS(), MatColoring
@*/
PetscErrorCode  ISColoringView(ISColoring iscoloring,PetscViewer viewer)
{
  PetscInt       i;
  PetscErrorCode ierr;
  PetscBool      iascii;
  IS             *is;

  PetscFunctionBegin;
  PetscValidPointer(iscoloring,1);
  if (!viewer) {
    ierr = PetscViewerASCIIGetStdout(iscoloring->comm,&viewer);CHKERRQ(ierr);
  }
  PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);

  ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
  if (iascii) {
    MPI_Comm    comm;
    PetscMPIInt size,rank;

    ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr);
    ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
    ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
    ierr = PetscViewerASCIIPrintf(viewer,"ISColoring Object: %d MPI processes\n",size);CHKERRQ(ierr);
    ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);CHKERRQ(ierr);
    ierr = PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of colors %d\n",rank,iscoloring->n);CHKERRQ(ierr);
    ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
    ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);CHKERRQ(ierr);
  } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for ISColoring",((PetscObject)viewer)->type_name);

  ierr = ISColoringGetIS(iscoloring,PETSC_IGNORE,&is);CHKERRQ(ierr);
  for (i=0; i<iscoloring->n; i++) {
    ierr = ISView(iscoloring->is[i],viewer);CHKERRQ(ierr);
  }
  ierr = ISColoringRestoreIS(iscoloring,&is);CHKERRQ(ierr);
  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);
}

/*@C
   MatGetColoring - Gets a coloring for a matrix, from its sparsity structure,
      to reduce the number of function evaluations needed to compute a sparse Jacobian via differencing.

   Collective on Mat

   Input Parameters:
.  mat - the matrix
.  type - type of coloring, one of the following:
$      MATCOLORINGNATURAL - natural (one color for each column, very slow)
$      MATCOLORINGSL - smallest-last
$      MATCOLORINGLF - largest-first
$      MATCOLORINGID - incidence-degree

   Output Parameters:
.   iscoloring - the coloring

   Options Database Keys:
   To specify the coloring through the options database, use one of
   the following
$    -mat_coloring_type natural, -mat_coloring_type sl, -mat_coloring_type lf,
$    -mat_coloring_type id
   To see the coloring use
$    -mat_coloring_view

   Level: intermediate

   Notes:
$    A graph coloring C(A)  is a division of vertices  so that two vertices of the same color do not share any common edges.
$    A suitable coloring for a  smoother  is simply C(A).
$    A suitable coloring for efficient Jacobian computation is a division of the columns so that two columns of the same color do not share any common rows.
$         This corresponds to C(A^{T} A).  This is what MatGetColoring() computes.

   The user can define additional colorings; see MatColoringRegisterDynamic().

   For parallel matrices currently converts to sequential matrix and uses the sequential coloring
   on that.

   The colorings SL, LF, and ID are obtained via the Minpack software that was
   converted to C using f2c.

   For BAIJ matrices this colors the blocks. The true number of colors would be block size times the number of colors
   returned here.

   References:
$     Thomas F. Coleman and Jorge J. More, Estimation of Sparse {J}acobian Matrices and Graph Coloring Problems,
$         SIAM Journal on Numerical Analysis, 1983, pages 187-209, volume 20
$     Jorge J. Mor\'{e} and Danny C. Sorenson and  Burton S. Garbow and Kenneth E. Hillstrom, The {MINPACK} Project,
$         Sources and Development of Mathematical Software, Wayne R. Cowell editor, 1984, pages 88-111

.keywords: matrix, get, coloring

.seealso:  MatGetColoringTypeFromOptions(), MatColoringRegisterDynamic(), MatFDColoringCreate(),
           SNESDefaultComputeJacobianColor()
@*/
PetscErrorCode  MatGetColoring(Mat mat,MatColoringType type,ISColoring *iscoloring)
{
  PetscBool      flag;
  PetscErrorCode ierr,(*r)(Mat,MatColoringType,ISColoring *);
  char           tname[PETSC_MAX_PATH_LEN];
  MPI_Comm       comm;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
  PetscValidPointer(iscoloring,3);
  if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
  if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

  /* look for type on command line */
  if (!MatColoringRegisterAllCalled) {ierr = MatColoringRegisterAll(PETSC_NULL);CHKERRQ(ierr);}
  ierr = PetscOptionsGetString(((PetscObject)mat)->prefix,"-mat_coloring_type",tname,256,&flag);CHKERRQ(ierr);
  if (flag) { type = tname; }

  ierr = PetscObjectGetComm((PetscObject)mat,&comm);CHKERRQ(ierr);
  ierr = PetscFunctionListFind(comm,MatColoringList,type,PETSC_TRUE,(void (**)(void)) &r);CHKERRQ(ierr);
  if (!r) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Unknown or unregistered type: %s",type);

  ierr = PetscLogEventBegin(MAT_GetColoring,mat,0,0,0);CHKERRQ(ierr);
  ierr = (*r)(mat,type,iscoloring);CHKERRQ(ierr);
  ierr = PetscLogEventEnd(MAT_GetColoring,mat,0,0,0);CHKERRQ(ierr);

  ierr = PetscInfo1(mat,"Number of colors %d\n",(*iscoloring)->n);CHKERRQ(ierr);
  flag = PETSC_FALSE;
  ierr = PetscOptionsGetBool(PETSC_NULL,"-mat_coloring_view",&flag,PETSC_NULL);CHKERRQ(ierr);
  if (flag) {
    PetscViewer viewer;
    ierr = PetscViewerASCIIGetStdout((*iscoloring)->comm,&viewer);CHKERRQ(ierr);
    ierr = ISColoringView(*iscoloring,viewer);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

/*@C
   KSPMonitorTrueResidualNorm - Prints the true residual norm as well as the preconditioned
   residual norm at each iteration of an iterative solver.

   Collective on KSP

   Input Parameters:
+  ksp   - iterative context
.  n     - iteration number
.  rnorm - 2-norm (preconditioned) residual value (may be estimated).
-  dummy - an ASCII PetscViewer

   Options Database Key:
.  -ksp_monitor_true_residual - Activates KSPMonitorTrueResidualNorm()

   Notes:
   When using right preconditioning, these values are equivalent.

   Level: intermediate

.keywords: KSP, default, monitor, residual

.seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(),KSPMonitorTrueResidualMaxNorm()
@*/
PetscErrorCode  KSPMonitorTrueResidualNorm(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
{
  PetscErrorCode ierr;
  Vec            resid;
  PetscReal      truenorm,bnorm;
  PetscViewer    viewer = (PetscViewer)dummy;
  char           normtype[256];

  PetscFunctionBegin;
  PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,4);
  ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);CHKERRQ(ierr);
  if (n == 0 && ((PetscObject)ksp)->prefix) {
    ierr = PetscViewerASCIIPrintf(viewer,"  Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);CHKERRQ(ierr);
  }
  ierr = KSPBuildResidual(ksp,NULL,NULL,&resid);CHKERRQ(ierr);
  ierr = VecNorm(resid,NORM_2,&truenorm);CHKERRQ(ierr);
  ierr = VecDestroy(&resid);CHKERRQ(ierr);
  ierr = VecNorm(ksp->vec_rhs,NORM_2,&bnorm);CHKERRQ(ierr);
  ierr = PetscStrncpy(normtype,KSPNormTypes[ksp->normtype],sizeof(normtype));CHKERRQ(ierr);
  ierr = PetscStrtolower(normtype);CHKERRQ(ierr);
  ierr = PetscViewerASCIIPrintf(viewer,"%3D KSP %s resid norm %14.12e true resid norm %14.12e ||r(i)||/||b|| %14.12e\n",n,normtype,(double)rnorm,(double)truenorm,(double)(truenorm/bnorm));CHKERRQ(ierr);
  ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);CHKERRQ(ierr);
  PetscFu