int DA::createActiveMatrix(Mat &M, MatType mtype, unsigned int dof) {
    // first determine the size ...
    unsigned int sz = 0;
    if(m_bIamActive) {
      sz = dof*(m_uiNodeSize + m_uiBoundaryNodeSize);

      // now create the PETSc Mat
      PetscBool isAij, isAijSeq, isAijPrl, isSuperLU, isSuperLU_Dist;
      PetscStrcmp(mtype,MATAIJ,&isAij);
      PetscStrcmp(mtype,MATSEQAIJ,&isAijSeq);
      PetscStrcmp(mtype,MATMPIAIJ,&isAijPrl);
      isSuperLU = PETSC_FALSE; //PetscStrcmp(mtype,MATSUPERLU,&isSuperLU);
      isSuperLU_Dist = PETSC_FALSE; //PetscStrcmp(mtype,MATSUPERLU_DIST,&isSuperLU_Dist);

      MatCreate(m_mpiCommActive, &M);
      MatSetSizes(M, sz,sz, PETSC_DECIDE, PETSC_DECIDE);
      MatSetType(M,mtype);

      if(isAij || isAijSeq || isAijPrl || isSuperLU || isSuperLU_Dist) {
        if(m_iNpesActive > 1) {
          MatMPIAIJSetPreallocation(M, 53*dof , PETSC_NULL, 53*dof , PETSC_NULL);
        }else {
          MatSeqAIJSetPreallocation(M, 53*dof , PETSC_NULL);
        }
      }
    }//end if active

    return 0;
  }//end function

PETSC_EXTERN PetscErrorCode TaoLMVMGetH0(Tao tao, Mat *H0)
{
  TAO_LMVM       *lmP;
  TAO_BLMVM      *blmP;
  const TaoType  type;
  PetscBool      is_lmvm, is_blmvm;
  Mat            M;

  PetscErrorCode ierr;

  ierr = TaoGetType(tao, &type);CHKERRQ(ierr);
  ierr = PetscStrcmp(type, TAOLMVM,  &is_lmvm);CHKERRQ(ierr);
  ierr = PetscStrcmp(type, TAOBLMVM, &is_blmvm);CHKERRQ(ierr);

  if (is_lmvm) {
    lmP = (TAO_LMVM *)tao->data;
    M = lmP->M;
  } else if (is_blmvm) {
    blmP = (TAO_BLMVM *)tao->data;
    M = blmP->M;
  } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONGSTATE, "This routine applies to TAO_LMVM and TAO_BLMVM.");

  ierr = MatLMVMGetH0(M, H0);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

int TaoLinearSolverPetsc::GetObjFcn(double *o_fcn)
{
  const KSPType ktype;
  int info;
  PetscTruth flg;

  PetscFunctionBegin;

  info = KSPGetType(ksp, &ktype); CHKERRQ(info);

  info = PetscStrcmp((char *)ktype, KSPNASH, &flg); CHKERRQ(info);
  if (flg == PETSC_TRUE) { 	
    info = KSPNASHGetObjFcn(ksp, o_fcn); CHKERRQ(info);
  }

  info = PetscStrcmp((char *)ktype, KSPSTCG, &flg); CHKERRQ(info);
  if (flg == PETSC_TRUE) { 	
    info = KSPSTCGGetObjFcn(ksp, o_fcn); CHKERRQ(info);
  }

  info = PetscStrcmp((char *)ktype, KSPGLTR, &flg); CHKERRQ(info);
  if (flg == PETSC_TRUE) { 	
    info = KSPGLTRGetObjFcn(ksp, o_fcn); CHKERRQ(info);
  }

  PetscFunctionReturn(0);
}

static PetscErrorCode PetscImageListAdd(const char filename[],const char ext[],PetscInt count)
{
  PetscErrorCode  ierr;
  PetscImageList  image,oimage = SAWs_images;
  PetscBool       flg;

  PetscFunctionBegin;
  if (oimage) {
    ierr = PetscStrcmp(filename,oimage->filename,&flg);CHKERRQ(ierr);
    if (flg) {
      oimage->count = count;
      PetscFunctionReturn(0);
    }
    while (oimage->next) {
      oimage = oimage->next;
      ierr = PetscStrcmp(filename,oimage->filename,&flg);CHKERRQ(ierr);
      if (flg) {
        oimage->count = count;
        PetscFunctionReturn(0);
      }
    }
    ierr = PetscNew(&image);CHKERRQ(ierr);
    oimage->next = image;
  } else {
    ierr = PetscRegisterFinalize(PetscImageListDestroy);CHKERRQ(ierr);
    ierr = PetscNew(&image);CHKERRQ(ierr);
    SAWs_images = image;
  }
  ierr = PetscStrallocpy(filename,&image->filename);CHKERRQ(ierr);
  ierr = PetscStrallocpy(ext,&image->ext);CHKERRQ(ierr);
  image->count = count;
  PetscFunctionReturn(0);
}

/*@C
   PetscMallocDump - Dumps the allocated memory blocks to a file. The information
   printed is: size of space (in bytes), address of space, id of space,
   file in which space was allocated, and line number at which it was
   allocated.

   Collective on PETSC_COMM_WORLD

   Input Parameter:
.  fp  - file pointer.  If fp is NULL, stdout is assumed.

   Options Database Key:
.  -malloc_dump - Dumps unfreed memory during call to PetscFinalize()

   Level: intermediate

   Fortran Note:
   The calling sequence in Fortran is PetscMallocDump(integer ierr)
   The fp defaults to stdout.

   Notes:
    uses MPI_COMM_WORLD, because this may be called in PetscFinalize() after PETSC_COMM_WORLD
          has been freed.

   Concepts: memory usage
   Concepts: memory bleeding
   Concepts: bleeding memory

.seealso:  PetscMallocGetCurrentUsage(), PetscMallocDumpLog()
@*/
PetscErrorCode  PetscMallocDump(FILE *fp)
{
  TRSPACE        *head;
  PetscInt       libAlloc = 0;
  PetscErrorCode ierr;
  PetscMPIInt    rank;

  PetscFunctionBegin;
  ierr = MPI_Comm_rank(MPI_COMM_WORLD,&rank);CHKERRQ(ierr);
  if (!fp) fp = PETSC_STDOUT;
  head = TRhead;
  while (head) {
    PetscBool isLib;

    ierr = PetscStrcmp(head->functionname, "PetscDLLibraryOpen", &isLib);CHKERRQ(ierr);
    libAlloc += head->size;
    head = head->next;
  }
  if (TRallocated - libAlloc > 0) fprintf(fp,"[%d]Total space allocated %.0f bytes\n",rank,(PetscLogDouble)TRallocated);
  head = TRhead;
  while (head) {
    PetscBool isLib;

    ierr = PetscStrcmp(head->functionname, "PetscDLLibraryOpen", &isLib);CHKERRQ(ierr);
    if (!isLib) {
      fprintf(fp,"[%2d]%.0f bytes %s() line %d in %s\n",rank,(PetscLogDouble)head->size,head->functionname,head->lineno,head->filename);
#if defined(PETSC_USE_DEBUG)
      ierr = PetscStackPrint(&head->stack,fp);CHKERRQ(ierr);
#endif
    }
    head = head->next;
  }
  PetscFunctionReturn(0);
}

/*@C
    PetscFOpen - Has the first process in the communicator open a file;
    all others do nothing.

    Logically Collective on MPI_Comm

    Input Parameters:
+   comm - the communicator
.   name - the filename
-   mode - the mode for fopen(), usually "w"

    Output Parameter:
.   fp - the file pointer

    Level: developer

    Notes:
       NULL (0), "stderr" or "stdout" may be passed in as the filename

    Fortran Note:
    This routine is not supported in Fortran.

    Concepts: opening ASCII file
    Concepts: files^opening ASCII

.seealso: PetscFClose(), PetscSynchronizedFGets(), PetscSynchronizedPrintf(), PetscSynchronizedFlush(),
          PetscFPrintf()
@*/
PetscErrorCode  PetscFOpen(MPI_Comm comm,const char name[],const char mode[],FILE **fp)
{
  PetscErrorCode ierr;
  PetscMPIInt    rank;
  FILE           *fd;
  char           fname[PETSC_MAX_PATH_LEN],tname[PETSC_MAX_PATH_LEN];

  PetscFunctionBegin;
  ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
  if (!rank) {
    PetscBool isstdout,isstderr;
    ierr = PetscStrcmp(name,"stdout",&isstdout);CHKERRQ(ierr);
    ierr = PetscStrcmp(name,"stderr",&isstderr);CHKERRQ(ierr);
    if (isstdout || !name) fd = PETSC_STDOUT;
    else if (isstderr) fd = PETSC_STDERR;
    else {
      PetscBool devnull;
      ierr = PetscStrreplace(PETSC_COMM_SELF,name,tname,PETSC_MAX_PATH_LEN);CHKERRQ(ierr);
      ierr = PetscFixFilename(tname,fname);CHKERRQ(ierr);
      ierr = PetscStrbeginswith(fname,"/dev/null",&devnull);CHKERRQ(ierr);
      if (devnull) {
        ierr = PetscStrcpy(fname,"/dev/null");CHKERRQ(ierr);
      }
      ierr = PetscInfo1(0,"Opening file %s\n",fname);CHKERRQ(ierr);
      fd   = fopen(fname,mode);
      if (!fd) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN,"Unable to open file %s\n",fname);
    }
  } else fd = 0;
  *fp = fd;
  PetscFunctionReturn(0);
}

EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "SNESMSSetType_MS"
PetscErrorCode  SNESMSSetType_MS(SNES snes,SNESMSType mstype)
{
  PetscErrorCode    ierr;
  SNES_MS           *ms = (SNES_MS*)snes->data;
  SNESMSTableauLink link;
  PetscBool         match;

  PetscFunctionBegin;
  if (ms->tableau) {
    ierr = PetscStrcmp(ms->tableau->name,mstype,&match);CHKERRQ(ierr);
    if (match) PetscFunctionReturn(0);
  }
  for (link = SNESMSTableauList; link; link=link->next) {
    ierr = PetscStrcmp(link->tab.name,mstype,&match);CHKERRQ(ierr);
    if (match) {
      ierr = SNESReset_MS(snes);CHKERRQ(ierr);
      ms->tableau = &link->tab;
      PetscFunctionReturn(0);
    }
  }
  SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",mstype);
  PetscFunctionReturn(0);
}

/*@C
  DMPlexRemoveLabel - Remove the label from this mesh

  Not Collective

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

  Output Parameter:
. label - The DMLabel, or NULL if the label is absent

  Level: developer

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

  PetscFunctionBegin;
  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
  ierr   = DMPlexHasLabel(dm, name, &hasLabel);CHKERRQ(ierr);
  *label = NULL;
  if (!hasLabel) PetscFunctionReturn(0);
  while (next) {
    ierr = PetscStrcmp(name, next->label->name, &hasLabel);CHKERRQ(ierr);
    if (hasLabel) {
      if (last) last->next   = next->next;
      else      mesh->labels = next->next;
      next->next = NULL;
      *label     = next->label;
      ierr = PetscFree(next);CHKERRQ(ierr);
      break;
    }
    last = next;
    next = next->next;
  }
  PetscFunctionReturn(0);
}

/*@
   MatAXPY - Computes Y = a*X + Y.

   Logically  Collective on Mat

   Input Parameters:
+  a - the scalar multiplier
.  X - the first matrix
.  Y - the second matrix
-  str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN
         or SUBSET_NONZERO_PATTERN (nonzeros of X is a subset of Y's)

   Level: intermediate

.keywords: matrix, add

.seealso: MatAYPX()
 @*/
PetscErrorCode MatAXPY(Mat Y,PetscScalar a,Mat X,MatStructure str)
{
  PetscErrorCode ierr;
  PetscInt       m1,m2,n1,n2;
  PetscBool      sametype;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(X,MAT_CLASSID,3);
  PetscValidHeaderSpecific(Y,MAT_CLASSID,1);
  PetscValidLogicalCollectiveScalar(Y,a,2);
  ierr = MatGetSize(X,&m1,&n1);CHKERRQ(ierr);
  ierr = MatGetSize(Y,&m2,&n2);CHKERRQ(ierr);
  if (m1 != m2 || n1 != n2) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Non conforming matrix add: %D %D %D %D",m1,m2,n1,n2);

  ierr = PetscStrcmp(((PetscObject)X)->type_name,((PetscObject)Y)->type_name,&sametype);CHKERRQ(ierr);
  ierr = PetscLogEventBegin(MAT_AXPY,Y,0,0,0);CHKERRQ(ierr);
  if (Y->ops->axpy && sametype) {
    ierr = (*Y->ops->axpy)(Y,a,X,str);CHKERRQ(ierr);
  } else {
    ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr);
  }
  ierr = PetscLogEventEnd(MAT_AXPY,Y,0,0,0);CHKERRQ(ierr);
#if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
  if (Y->valid_GPU_matrix != PETSC_OFFLOAD_UNALLOCATED) {
    Y->valid_GPU_matrix = PETSC_OFFLOAD_CPU;
  }
#endif
  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);
}

/*@C
   DMGetNamedLocalVector - get access to a named, persistent local vector

   Not Collective

   Input Arguments:
+  dm - DM to hold named vectors
-  name - unique name for Vec

   Output Arguments:
.  X - named Vec

   Level: developer

   Note: If a Vec with the given name does not exist, it is created.

.seealso: DMGetNamedGlobalVector(),DMRestoreNamedLocalVector()
@*/
PetscErrorCode DMGetNamedLocalVector(DM dm,const char *name,Vec *X)
{
  PetscErrorCode ierr;
  DMNamedVecLink link;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
  PetscValidCharPointer(name,2);
  PetscValidPointer(X,3);
  for (link=dm->namedlocal; link; link=link->next) {
    PetscBool match;
    ierr = PetscStrcmp(name,link->name,&match);CHKERRQ(ierr);
    if (match) {
      if (link->status != DMVEC_STATUS_IN) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Vec name '%s' already checked out",name);
      goto found;
    }
  }

  /* Create the Vec */
  ierr           = PetscNew(&link);CHKERRQ(ierr);
  ierr           = PetscStrallocpy(name,&link->name);CHKERRQ(ierr);
  ierr           = DMCreateLocalVector(dm,&link->X);CHKERRQ(ierr);
  link->next     = dm->namedlocal;
  dm->namedlocal = link;

found:
  *X           = link->X;
  link->status = DMVEC_STATUS_OUT;
  PetscFunctionReturn(0);
}

PetscErrorCode PetscViewerHDF5OpenGroup(PetscViewer viewer, hid_t *fileId, hid_t *groupId)
{
  hid_t          file_id, group;
  const char     *groupName = NULL;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = PetscViewerHDF5GetFileId(viewer, &file_id);CHKERRQ(ierr);
  ierr = PetscViewerHDF5GetGroup(viewer, &groupName);CHKERRQ(ierr);
  /* Open group */
  if (groupName) {
    PetscBool root;

    ierr = PetscStrcmp(groupName, "/", &root);CHKERRQ(ierr);
    if (!root && !H5Lexists(file_id, groupName, H5P_DEFAULT)) {
#if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800)
      group = H5Gcreate2(file_id, groupName, 0, H5P_DEFAULT, H5P_DEFAULT);
#else /* deprecated HDF5 1.6 API */
      group = H5Gcreate(file_id, groupName, 0);
#endif
      if (group < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "Could not create group %s", groupName);
      ierr = H5Gclose(group);CHKERRQ(ierr);
    }
#if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800)
    group = H5Gopen2(file_id, groupName, H5P_DEFAULT);
#else
    group = H5Gopen(file_id, groupName);
#endif
    if (group < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "Could not open group %s", groupName);
  } else group = file_id;

  *fileId  = file_id;
  *groupId = group;
  PetscFunctionReturn(0);
}

static PetscErrorCode DMCreateMatrix_Shell(DM dm,Mat *J)
{
    PetscErrorCode ierr;
    DM_Shell       *shell = (DM_Shell*)dm->data;
    Mat            A;

    PetscFunctionBegin;
    PetscValidHeaderSpecific(dm,DM_CLASSID,1);
    PetscValidPointer(J,3);
    if (!shell->A) {
        if (shell->Xglobal) {
            PetscInt m,M;
            ierr = PetscInfo(dm,"Naively creating matrix using global vector distribution without preallocation\n");
            CHKERRQ(ierr);
            ierr = VecGetSize(shell->Xglobal,&M);
            CHKERRQ(ierr);
            ierr = VecGetLocalSize(shell->Xglobal,&m);
            CHKERRQ(ierr);
            ierr = MatCreate(PetscObjectComm((PetscObject)dm),&shell->A);
            CHKERRQ(ierr);
            ierr = MatSetSizes(shell->A,m,m,M,M);
            CHKERRQ(ierr);
            ierr = MatSetType(shell->A,dm->mattype);
            CHKERRQ(ierr);
            ierr = MatSetUp(shell->A);
            CHKERRQ(ierr);
        } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMShellSetMatrix(), DMShellSetCreateMatrix(), or provide a vector");
    }
    A = shell->A;
    /* the check below is tacky and incomplete */
    if (dm->mattype) {
        PetscBool flg,aij,seqaij,mpiaij;
        ierr = PetscObjectTypeCompare((PetscObject)A,dm->mattype,&flg);
        CHKERRQ(ierr);
        ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&seqaij);
        CHKERRQ(ierr);
        ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&mpiaij);
        CHKERRQ(ierr);
        ierr = PetscStrcmp(dm->mattype,MATAIJ,&aij);
        CHKERRQ(ierr);
        if (!flg) {
            if (!(aij && (seqaij || mpiaij))) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_NOTSAMETYPE,"Requested matrix of type %s, but only %s available",dm->mattype,((PetscObject)A)->type_name);
        }
    }
    if (((PetscObject)A)->refct < 2) { /* We have an exclusive reference so we can give it out */
        ierr = PetscObjectReference((PetscObject)A);
        CHKERRQ(ierr);
        ierr = MatZeroEntries(A);
        CHKERRQ(ierr);
        *J   = A;
    } else {                      /* Need to create a copy, could use MAT_SHARE_NONZERO_PATTERN in most cases */
        ierr = MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,J);
        CHKERRQ(ierr);
        ierr = MatZeroEntries(*J);
        CHKERRQ(ierr);
    }
    PetscFunctionReturn(0);
}

/*@C
   PetscDrawSetType - Builds graphics object for a particular implementation

   Collective on PetscDraw

   Input Parameter:
+  draw      - the graphics context
-  type      - for example, PETSC_DRAW_X

   Options Database Command:
.  -draw_type  <type> - Sets the type; use -help for a list
    of available methods (for instance, x)

   Level: intermediate

   Notes:
   See "petsc/include/petscdraw.h" for available methods (for instance,
   PETSC_DRAW_X)

   Concepts: drawing^X windows
   Concepts: X windows^graphics
   Concepts: drawing^Microsoft Windows

.seealso: PetscDrawSetFromOptions(), PetscDrawCreate(), PetscDrawDestroy()
@*/
PetscErrorCode  PetscDrawSetType(PetscDraw draw,PetscDrawType type)
{
    PetscErrorCode ierr,(*r)(PetscDraw);
    PetscBool      match;
    PetscBool      flg=PETSC_FALSE;

    PetscFunctionBegin;
    PetscValidHeaderSpecific(draw,PETSC_DRAW_CLASSID,1);
    PetscValidCharPointer(type,2);

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

    /*  User requests no graphics */
    ierr = PetscOptionsHasName(NULL,"-nox",&flg);
    CHKERRQ(ierr);

    /*
       This is not ideal, but it allows codes to continue to run if X graphics
     was requested but is not installed on this machine. Mostly this is for
     testing.
     */
#if !defined(PETSC_HAVE_X)
    if (!flg) {
        ierr = PetscStrcmp(type,PETSC_DRAW_X,&match);
        CHKERRQ(ierr);
        if (match) {
            PetscBool dontwarn = PETSC_TRUE;
            flg  = PETSC_TRUE;
            ierr = PetscOptionsHasName(NULL,"-nox_warning",&dontwarn);
            CHKERRQ(ierr);
            if (!dontwarn) (*PetscErrorPrintf)("PETSc installed without X windows on this machine\nproceeding without graphics\n");
        }
    }
#endif
    if (flg) type = PETSC_DRAW_NULL;

    if (draw->data) {
        /* destroy the old private PetscDraw context */
        ierr               = (*draw->ops->destroy)(draw);
        CHKERRQ(ierr);
        draw->ops->destroy = NULL;
        draw->data         = 0;
    }

    ierr =  PetscFunctionListFind(PetscDrawList,type,&r);
    CHKERRQ(ierr);
    if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown PetscDraw type given: %s",type);
    ierr       = PetscObjectChangeTypeName((PetscObject)draw,type);
    CHKERRQ(ierr);
    draw->data = 0;
    ierr       = (*r)(draw);
    CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

/*@C
    PetscPushJSONValue -  Puts a "key" : "value" pair onto a string

    Input Parameters:
+   buffer - the char array where the value will be put
.   key - the key value to be set
.   value - the value associated with the key
-   bufflen - the size of the buffer (currently ignored)

    Level: advanced

    Notes:
    Ignores lengths so can cause buffer overflow
@*/
PetscErrorCode PetscPushJSONValue(char buff[],const char key[],const char value[],size_t bufflen)
{
  PetscErrorCode ierr;
  size_t         len;
  PetscBool      special;

  PetscFunctionBegin;
  ierr = PetscStrcmp(value,"null",&special);CHKERRQ(ierr);
  if (!special) {
    ierr = PetscStrcmp(value,"true",&special);CHKERRQ(ierr);
  }
  if (!special) {
    ierr = PetscStrcmp(value,"false",&special);CHKERRQ(ierr);
  }
  if (!special) {
    PetscInt i;

    ierr    = PetscStrlen(value,&len);CHKERRQ(ierr);
    special = PETSC_TRUE;
    for (i=0; i<(int)len; i++) {
      if (!isdigit(value[i])) {
        special = PETSC_FALSE;
        break;
      }
    }
  }

  ierr = PetscStrcat(buff,"\"");CHKERRQ(ierr);
  ierr = PetscStrcat(buff,key);CHKERRQ(ierr);
  ierr = PetscStrcat(buff,"\":");CHKERRQ(ierr);
  if (!special) {
    ierr = PetscStrcat(buff,"\"");CHKERRQ(ierr);
  }
  ierr = PetscStrcat(buff,value);CHKERRQ(ierr);
  if (!special) {
    ierr = PetscStrcat(buff,"\"");CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

dErr dUnitsFindUnit(dUnits un,const char *name,dUnit *unit)
{
  dErr err;
  dFunctionBegin;
  *unit = NULL;
  for (dInt i=0; i<un->nalloc; i++) {
    dBool flg;
    dUnit t = un->list[i];
    err = PetscStrcmp(t->longname,name,&flg);dCHK(err);
    if (flg) {*unit = t; break;}
  }
  dFunctionReturn(0);
}

PetscErrorCode  TSRKSetType_RK(TS ts,TSRKType rktype)
{
  TS_RK     *rk = (TS_RK*)ts->data;
  PetscErrorCode ierr;
  PetscBool      match;
  RKTableauLink link;

  PetscFunctionBegin;
  if (rk->tableau) {
    ierr = PetscStrcmp(rk->tableau->name,rktype,&match);CHKERRQ(ierr);
    if (match) PetscFunctionReturn(0);
  }
  for (link = RKTableauList; link; link=link->next) {
    ierr = PetscStrcmp(link->tab.name,rktype,&match);CHKERRQ(ierr);
    if (match) {
      ierr = TSReset_RK(ts);CHKERRQ(ierr);
      rk->tableau = &link->tab;
      PetscFunctionReturn(0);
    }
  }
  SETERRQ1(PetscObjectComm((PetscObject)ts),PETSC_ERR_ARG_UNKNOWN_TYPE,"Could not find '%s'",rktype);
  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);
}

/*MC
   PetscFunctionListAdd - Given a routine and a string id, saves that routine in the
   specified registry.

   Synopsis:
   #include <petscsys.h>
   PetscErrorCode PetscFunctionListAdd(PetscFunctionList *flist,const char name[],void (*fptr)(void))

   Not Collective

   Input Parameters:
+  flist - pointer to function list object
.  name - string to identify routine
-  fptr - function pointer

   Notes:
   To remove a registered routine, pass in a NULL fptr.

   Users who wish to register new classes for use by a particular PETSc
   component (e.g., SNES) should generally call the registration routine
   for that particular component (e.g., SNESRegister()) instead of
   calling PetscFunctionListAdd() directly.

    Level: developer

.seealso: PetscFunctionListDestroy(), SNESRegister(), KSPRegister(),
          PCRegister(), TSRegister(), PetscFunctionList, PetscObjectComposeFunction()
M*/
PETSC_EXTERN PetscErrorCode PetscFunctionListAdd_Private(PetscFunctionList *fl,const char name[],void (*fnc)(void))
{
  PetscFunctionList entry,ne;
  PetscErrorCode    ierr;

  PetscFunctionBegin;
  if (!*fl) {
    ierr           = PetscNew(&entry);CHKERRQ(ierr);
    ierr           = PetscStrallocpy(name,&entry->name);CHKERRQ(ierr);
    entry->routine = fnc;
    entry->next    = 0;
    *fl            = entry;

#if defined(PETSC_USE_DEBUG)
    /* add this new list to list of all lists */
    if (!dlallhead) {
      dlallhead        = *fl;
      (*fl)->next_list = 0;
    } else {
      ne               = dlallhead;
      dlallhead        = *fl;
      (*fl)->next_list = ne;
    }
#endif

  } else {
    /* search list to see if it is already there */
    ne = *fl;
    while (ne) {
      PetscBool founddup;

      ierr = PetscStrcmp(ne->name,name,&founddup);CHKERRQ(ierr);
      if (founddup) { /* found duplicate */
        ne->routine = fnc;
        PetscFunctionReturn(0);
      }
      if (ne->next) ne = ne->next;
      else break;
    }
    /* create new entry and add to end of list */
    ierr           = PetscNew(&entry);CHKERRQ(ierr);
    ierr           = PetscStrallocpy(name,&entry->name);CHKERRQ(ierr);
    entry->routine = fnc;
    entry->next    = 0;
    ne->next       = entry;
  }
  PetscFunctionReturn(0);
}

/*@C
   PetscObjectTypeCompare - Determines whether a PETSc object is of a particular type.

   Not Collective

   Input Parameters:
+  obj - any PETSc object, for example a Vec, Mat or KSP.
         This must be cast with a (PetscObject), for example,
         PetscObjectTypeCompare((PetscObject)mat);
-  type_name - string containing a type name

   Output Parameter:
.  same - PETSC_TRUE if they are the same, else PETSC_FALSE

   Level: intermediate

.seealso: VecGetType(), KSPGetType(), PCGetType(), SNESGetType()

   Concepts: comparing^object types
   Concepts: types^comparing
   Concepts: object type^comparpeing

@*/
PetscErrorCode  PetscObjectTypeCompare(PetscObject obj,const char type_name[],PetscBool  *same)
{
  PetscErrorCode ierr;

  PetscFunctionBegin;
  if (!obj) *same = PETSC_FALSE;
  else if (!type_name && !obj->type_name) *same = PETSC_TRUE;
  else if (!type_name || !obj->type_name) *same = PETSC_FALSE;
  else {
    PetscValidHeader(obj,1);
    PetscValidCharPointer(type_name,2);
    PetscValidPointer(same,3);
    ierr = PetscStrcmp((char*)(obj->type_name),type_name,same);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}