/*@C
   PetscDrawSetSave - Saves images produced in a PetscDraw into a file as a Gif file using AfterImage

   Collective on PetscDraw

   Input Parameter:
+  draw      - the graphics context
.  filename  - name of the file, if NULL uses name of draw object
-  movie - produce a movie of all the images

   Options Database Command:
+  -draw_save  <filename>
-  -draw_save_movie

   Level: intermediate

   Concepts: X windows^graphics

   Notes: You should call this BEFORE calling PetscDrawClear() and creating your image.

   Requires that PETSc be configured with the option --with-afterimage to save the images and ffmpeg must be in your path to make the movie

   If X windows generates an error message about X_CreateWindow() failing then Afterimage was installed without X windows. Reinstall Afterimage using the
   ./configure flags --x-includes=/pathtoXincludes --x-libraries=/pathtoXlibraries   For example under Mac OS X Mountain Lion --x-includes=/opt/X11/include -x-libraries=/opt/X11/lib


.seealso: PetscDrawSetFromOptions(), PetscDrawCreate(), PetscDrawDestroy()
@*/
PetscErrorCode  PetscDrawSetSave(PetscDraw draw,const char *filename,PetscBool movie)
{
    PetscErrorCode ierr;

    PetscFunctionBegin;
    PetscValidHeaderSpecific(draw,PETSC_DRAW_CLASSID,1);
    ierr = PetscFree(draw->savefilename);
    CHKERRQ(ierr);

    draw->savefilemovie = movie;
    if (filename && filename[0]) {
        ierr = PetscStrallocpy(filename,&draw->savefilename);
        CHKERRQ(ierr);
    } else {
        const char *name;
        ierr = PetscObjectGetName((PetscObject)draw,&name);
        CHKERRQ(ierr);
        ierr = PetscStrallocpy(name,&draw->savefilename);
        CHKERRQ(ierr);
    }
    if (draw->ops->setsave) {
        ierr = (*draw->ops->setsave)(draw,filename);
        CHKERRQ(ierr);
    }
    PetscFunctionReturn(0);
}

/*@C
  DMPlexAddBoundary - Add a boundary condition to the model

  Input Parameters:
+ dm          - The mesh object
. isEssential - Flag for an essential (Dirichlet) condition, as opposed to a natural (Neumann) condition
. name        - The BC name
. labelname   - The label defining constrained points
. field       - The field to constrain
. numcomps    - The number of constrained field components
. comps       - An array of constrained component numbers
. bcFunc      - A pointwise function giving boundary values
. numids      - The number of DMLabel ids for constrained points
. ids         - An array of ids for constrained points
- ctx         - An optional user context for bcFunc

  Options Database Keys:
+ -bc_<boundary name> <num> - Overrides the boundary ids
- -bc_<boundary name>_comp <num> - Overrides the boundary components

  Level: developer

.seealso: DMPlexGetBoundary()
@*/
PetscErrorCode DMPlexAddBoundary(DM dm, PetscBool isEssential, const char name[], const char labelname[], PetscInt field, PetscInt numcomps, const PetscInt *comps, void (*bcFunc)(), PetscInt numids, const PetscInt *ids, void *ctx)
{
  DM_Plex       *mesh = (DM_Plex *) dm->data;
  DMBoundary     b;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
  ierr = PetscNew(&b);CHKERRQ(ierr);
  ierr = PetscStrallocpy(name, (char **) &b->name);CHKERRQ(ierr);
  ierr = PetscStrallocpy(labelname, (char **) &b->labelname);CHKERRQ(ierr);
  ierr = PetscMalloc1(numcomps, &b->comps);CHKERRQ(ierr);
  if (numcomps) {ierr = PetscMemcpy(b->comps, comps, numcomps*sizeof(PetscInt));CHKERRQ(ierr);}
  ierr = PetscMalloc1(numids, &b->ids);CHKERRQ(ierr);
  if (numids) {ierr = PetscMemcpy(b->ids, ids, numids*sizeof(PetscInt));CHKERRQ(ierr);}
  if (b->labelname) {
    ierr = DMPlexGetLabel(dm, b->labelname, &b->label);CHKERRQ(ierr);
    if (!b->label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Label %s does not exist in this DM", b->labelname);
  }
  b->essential   = isEssential;
  b->field       = field;
  b->numcomps    = numcomps;
  b->func        = bcFunc;
  b->numids      = numids;
  b->ctx         = ctx;
  b->next        = mesh->boundary;
  mesh->boundary = b;
  PetscFunctionReturn(0);
}

PETSC_EXTERN PetscErrorCode PetscDrawCreate_TikZ(PetscDraw draw)
{
  PetscDraw_TikZ *win;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = PetscMemcpy(draw->ops,&DvOps,sizeof(DvOps));CHKERRQ(ierr);
  ierr = PetscNew(PetscDraw_TikZ,&win);CHKERRQ(ierr);
  ierr = PetscLogObjectMemory(draw,sizeof(PetscDraw_TikZ));CHKERRQ(ierr);

  draw->data = (void*) win;

  if (draw->title) {
    ierr = PetscStrallocpy(draw->title,&win->filename);CHKERRQ(ierr);
  } else {
    const char *fname;
    ierr = PetscObjectGetName((PetscObject)draw,&fname);CHKERRQ(ierr);
    ierr = PetscStrallocpy(fname,&win->filename);CHKERRQ(ierr);
  }
  ierr = PetscFOpen(PetscObjectComm((PetscObject)draw),win->filename,"w",&win->fd);CHKERRQ(ierr);
  ierr = PetscFPrintf(PetscObjectComm((PetscObject)draw),win->fd,TikZ_BEGIN_DOCUMENT);CHKERRQ(ierr);
  ierr = PetscFPrintf(PetscObjectComm((PetscObject)draw),win->fd,TikZ_BEGIN_FRAME);CHKERRQ(ierr);

  win->written = PETSC_FALSE;
  PetscFunctionReturn(0);
}

static PetscErrorCode BoundaryDuplicate(DMBoundary bd, DMBoundary *boundary)
{
  DMBoundary     b = bd, b2, bold = NULL;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  *boundary = NULL;
  for (; b; b = b->next, bold = b2) {
    ierr = PetscNew(&b2);CHKERRQ(ierr);
    ierr = PetscStrallocpy(b->name, (char **) &b2->name);CHKERRQ(ierr);
    ierr = PetscStrallocpy(b->labelname, (char **) &b2->labelname);CHKERRQ(ierr);
    ierr = PetscMalloc1(b->numids, &b2->ids);CHKERRQ(ierr);
    ierr = PetscMemcpy(b2->ids, b->ids, b->numids*sizeof(PetscInt));CHKERRQ(ierr);
    ierr = PetscMalloc1(b->numcomps, &b2->comps);CHKERRQ(ierr);
    ierr = PetscMemcpy(b2->comps, b->comps, b->numcomps*sizeof(PetscInt));CHKERRQ(ierr);
    b2->label     = NULL;
    b2->essential = b->essential;
    b2->field     = b->field;
    b2->numcomps  = b->numcomps;
    b2->func      = b->func;
    b2->numids    = b->numids;
    b2->ctx       = b->ctx;
    b2->next      = NULL;
    if (!*boundary) *boundary   = b2;
    if (bold)        bold->next = b2;
  }
  PetscFunctionReturn(0);
}

PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_matlab(Mat A,MatFactorType ftype,Mat *F)
{
    PetscErrorCode ierr;

    PetscFunctionBegin;
    if (A->cmap->N != A->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"matrix must be square");
    ierr                         = MatCreate(PetscObjectComm((PetscObject)A),F);
    CHKERRQ(ierr);
    ierr                         = MatSetSizes(*F,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
    CHKERRQ(ierr);
    ierr                         = PetscStrallocpy("matlab",&((PetscObject)*F)->type_name);
    CHKERRQ(ierr);
    ierr                         = MatSetUp(*F);
    CHKERRQ(ierr);

    (*F)->ops->destroy           = MatDestroy_matlab;
    (*F)->ops->getinfo           = MatGetInfo_External;
    (*F)->ops->lufactorsymbolic  = MatLUFactorSymbolic_Matlab;
    (*F)->ops->ilufactorsymbolic = MatLUFactorSymbolic_Matlab;

    ierr = PetscObjectComposeFunction((PetscObject)(*F),"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_seqaij_matlab);
    CHKERRQ(ierr);

    (*F)->factortype = ftype;
    ierr = PetscFree((*F)->solvertype);
    CHKERRQ(ierr);
    ierr = PetscStrallocpy(MATSOLVERMATLAB,&(*F)->solvertype);
    CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

static dErr dUnitsAssignName(dUnits un,const char *(*namer)(dUnit),const char *proposed,dInt n,const dReal expon[],char **assigned)
{
  dErr err;
  char buf[1024],*p = buf;
  dInt left = 1024;
  dFunctionBegin;
  if (proposed) {
    err = PetscStrallocpy(proposed,assigned);dCHK(err);
    dFunctionReturn(0);
  }
  buf[0] = 0;
  for (dInt i=0; i<n; i++) {
    const char *s;
    dUnit base;
    dInt len;
    if (expon[i] == 0) continue;
    err = dUnitsGetBase(un,i,&base);dCHK(err);
    s = namer(base);
    if (expon[i] == 1)
      len = snprintf(p,left,"%s ",s);
    else if (round(expon[i]) == expon[i])
      len = snprintf(p,left,"%s^%1.0f ",s,expon[i]);
    else
      len = snprintf(p,left,"%s^%f ",s,expon[i]);
    left -= len;
    p += len;
  }
  p[-1] = 0; // Kill trailing space
  err = PetscStrallocpy(buf,assigned);dCHK(err);
  dFunctionReturn(0);
}

PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_essl(Mat A,MatFactorType ftype,Mat *F)
{
  Mat            B;
  PetscErrorCode ierr;
  Mat_Essl       *essl;

  PetscFunctionBegin;
  if (A->cmap->N != A->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"matrix must be square");
  ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr);
  ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,A->rmap->n,A->cmap->n);CHKERRQ(ierr);
  ierr = PetscStrallocpy("essl",&((PetscObject)B)->type_name);CHKERRQ(ierr);
  ierr = MatSetUp(B);CHKERRQ(ierr);

  ierr = PetscNewLog(B,&essl);CHKERRQ(ierr);

  B->data                  = essl;
  B->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl;
  B->ops->destroy          = MatDestroy_Essl;
  B->ops->getinfo          = MatGetInfo_External;

  ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_essl);CHKERRQ(ierr);

  B->factortype = MAT_FACTOR_LU;
  ierr = PetscFree(B->solvertype);CHKERRQ(ierr);
  ierr = PetscStrallocpy(MATSOLVERESSL,&B->solvertype);CHKERRQ(ierr);

  *F            = B;
  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);
}

PetscErrorCode  PetscViewerDrawSetInfo(PetscViewer v,const char display[],const char title[],int x,int y,int w,int h)
{
  PetscErrorCode   ierr;
  PetscViewer_Draw *vdraw = (PetscViewer_Draw*)v->data;

  PetscFunctionBegin;
  vdraw->h = h;
  vdraw->w = w;
  ierr     = PetscStrallocpy(display,&vdraw->display);CHKERRQ(ierr);
  ierr     = PetscStrallocpy(title,&vdraw->title);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*@C
   PetscDrawCreate - Creates a graphics context.

   Collective on MPI_Comm

   Input Parameter:
+  comm - MPI communicator
.  display - X display when using X windows
.  title - optional title added to top of window
.  x,y - coordinates of lower left corner of window or PETSC_DECIDE
-  w, h - width and height of window or PETSC_DECIDE or PETSC_DRAW_HALF_SIZE, PETSC_DRAW_FULL_SIZE,
          or PETSC_DRAW_THIRD_SIZE or PETSC_DRAW_QUARTER_SIZE

   Output Parameter:
.  draw - location to put the PetscDraw context

   Level: beginner

   Concepts: graphics^creating context
   Concepts: drawing^creating context

.seealso: PetscDrawSetFromOptions(), PetscDrawDestroy(), PetscDrawSetType()
@*/
PetscErrorCode  PetscDrawCreate(MPI_Comm comm,const char display[],const char title[],int x,int y,int w,int h,PetscDraw *indraw)
{
    PetscDraw      draw;
    PetscErrorCode ierr;
    PetscReal      dpause;
    PetscBool      flag;

    PetscFunctionBegin;
    ierr = PetscDrawInitializePackage();
    CHKERRQ(ierr);
    *indraw = 0;
    ierr = PetscHeaderCreate(draw,_p_PetscDraw,struct _PetscDrawOps,PETSC_DRAW_CLASSID,"Draw","Graphics","Draw",comm,PetscDrawDestroy,0);
    CHKERRQ(ierr);

    draw->data    = 0;
    ierr          = PetscStrallocpy(title,&draw->title);
    CHKERRQ(ierr);
    ierr          = PetscStrallocpy(display,&draw->display);
    CHKERRQ(ierr);
    draw->x       = x;
    draw->y       = y;
    draw->w       = w;
    draw->h       = h;
    draw->pause   = 0.0;
    draw->coor_xl = 0.0;
    draw->coor_xr = 1.0;
    draw->coor_yl = 0.0;
    draw->coor_yr = 1.0;
    draw->port_xl = 0.0;
    draw->port_xr = 1.0;
    draw->port_yl = 0.0;
    draw->port_yr = 1.0;
    draw->popup   = 0;

    ierr = PetscOptionsGetReal(NULL,"-draw_pause",&dpause,&flag);
    CHKERRQ(ierr);
    if (flag) draw->pause = dpause;
    draw->savefilename  = NULL;
    draw->savefilemovie = PETSC_FALSE;
    draw->savefilecount = -1;

    ierr = PetscDrawSetCurrentPoint(draw,.5,.9);
    CHKERRQ(ierr);

    draw->boundbox_xl  = .5;
    draw->boundbox_xr  = .5;
    draw->boundbox_yl  = .9;
    draw->boundbox_yr  = .9;

    *indraw = draw;
    PetscFunctionReturn(0);
}

/*@C
    PetscDrawAxisSetLabels -  Sets the x and y axis labels.

    Not Collective (ignored on all processors except processor 0 of PetscDrawAxis)

    Input Parameters:
+   axis - the axis
.   top - the label at the top of the image
-   xlabel,ylabel - the labes for the x and y axis

    Notes: Must be called before PetscDrawAxisDraw() or PetscDrawLGDraw()
           There should be no newlines in the arguments

    Level: advanced

@*/
PetscErrorCode  PetscDrawAxisSetLabels(PetscDrawAxis axis,const char top[],const char xlabel[],const char ylabel[])
{
  PetscErrorCode ierr;

  PetscFunctionBegin;
  if (!axis) PetscFunctionReturn(0);
  ierr = PetscFree(axis->xlabel);CHKERRQ(ierr);
  ierr = PetscFree(axis->ylabel);CHKERRQ(ierr);
  ierr = PetscFree(axis->toplabel);CHKERRQ(ierr);
  ierr = PetscStrallocpy(xlabel,&axis->xlabel);CHKERRQ(ierr);
  ierr = PetscStrallocpy(ylabel,&axis->ylabel);CHKERRQ(ierr);
  ierr = PetscStrallocpy(top,&axis->toplabel);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PETSC_EXTERN PetscErrorCode PCCreate_ICC(PC pc)
{
  PetscErrorCode ierr;
  PC_ICC         *icc;

  PetscFunctionBegin;
  ierr = PetscNewLog(pc,&icc);CHKERRQ(ierr);

  ((PC_Factor*)icc)->fact = 0;

  ierr = PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)icc)->ordering);CHKERRQ(ierr);
  ierr = PetscStrallocpy(MATSOLVERPETSC,&((PC_Factor*)icc)->solvertype);CHKERRQ(ierr);
  ierr = MatFactorInfoInitialize(&((PC_Factor*)icc)->info);CHKERRQ(ierr);

  ((PC_Factor*)icc)->factortype  = MAT_FACTOR_ICC;
  ((PC_Factor*)icc)->info.levels = 0.;
  ((PC_Factor*)icc)->info.fill   = 1.0;
  icc->implctx                   = 0;

  ((PC_Factor*)icc)->info.dtcol       = PETSC_DEFAULT;
  ((PC_Factor*)icc)->info.shifttype   = (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE;
  ((PC_Factor*)icc)->info.shiftamount = 100.0*PETSC_MACHINE_EPSILON;
  ((PC_Factor*)icc)->info.zeropivot   = 100.0*PETSC_MACHINE_EPSILON;

  pc->data                     = (void*)icc;
  pc->ops->apply               = PCApply_ICC;
  pc->ops->applytranspose      = PCApply_ICC;
  pc->ops->setup               = PCSetup_ICC;
  pc->ops->reset               = PCReset_ICC;
  pc->ops->destroy             = PCDestroy_ICC;
  pc->ops->setfromoptions      = PCSetFromOptions_ICC;
  pc->ops->view                = PCView_ICC;
  pc->ops->getfactoredmatrix   = PCFactorGetMatrix_Factor;
  pc->ops->applysymmetricleft  = PCApplySymmetricLeft_ICC;
  pc->ops->applysymmetricright = PCApplySymmetricRight_ICC;

  ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C",PCFactorSetUpMatSolverPackage_Factor);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetMatSolverPackage_C",PCFactorGetMatSolverPackage_Factor);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetZeroPivot_C",PCFactorSetZeroPivot_Factor);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftType_C",PCFactorSetShiftType_Factor);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftAmount_C",PCFactorSetShiftAmount_Factor);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetLevels_C",PCFactorSetLevels_Factor);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetLevels_C",PCFactorGetLevels_Factor);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetFill_C",PCFactorSetFill_Factor);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatOrderingType_C",PCFactorSetMatOrderingType_Factor);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatSolverPackage_C",PCFactorSetMatSolverPackage_Factor);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetDropTolerance_C",PCFactorSetDropTolerance_ILU);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetUseInPlace_C",PCFactorGetUseInPlace_ICC);CHKERRQ(ierr);
  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
   PetscDrawSetSave - Saves images produced in a PetscDraw into a file

   Collective on PetscDraw

   Input Parameter:
+  draw      - the graphics context
.  filename  - name of the file, if .ext then uses name of draw object plus .ext using .ext to determine the image type
-  movieext  - if not NULL, produces a movie of all the images

   Options Database Command:
+  -draw_save  <filename>  - filename could be name.ext or .ext (where .ext determines the type of graphics file to save, for example .png)
.  -draw_save_movie <.ext> - saves a movie to filename.ext
.  -draw_save_final_image [optional filename] - saves the final image displayed in a window
-  -draw_save_single_file - saves each new image in the same file, normally each new image is saved in a new file with filename/filename_%d.ext

   Level: intermediate

   Concepts: X windows^graphics

   Notes: You should call this BEFORE creating your image and calling PetscDrawSave().
   The supported image types are .png, .gif, .jpg, and .ppm (PETSc chooses the default in that order).
   Support for .png images requires configure --with-libpng.
   Support for .gif images requires configure --with-giflib.
   Support for .jpg images requires configure --with-libjpeg.
   Support for .ppm images is built-in. The PPM format has no compression (640x480 pixels ~ 900 KiB).
   The ffmpeg utility must be in your path to make the movie.

.seealso: PetscDrawSetFromOptions(), PetscDrawCreate(), PetscDrawDestroy(), PetscDrawSetSaveFinalImage()
@*/
PetscErrorCode  PetscDrawSetSave(PetscDraw draw,const char filename[],const char movieext[])
{
  const char     *savename = NULL;
  const char     *imageext = NULL;
  char           buf[PETSC_MAX_PATH_LEN];
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(draw,PETSC_DRAW_CLASSID,1);
  if (filename) PetscValidCharPointer(filename,2);
  if (movieext) PetscValidCharPointer(movieext,2);

  /* determine save filename and image extension */
  if (filename && filename[0]) {
    ierr = PetscStrchr(filename,'.',(char **)&imageext);CHKERRQ(ierr);
    if (!imageext) savename = filename;
    else if (imageext != filename) {
      size_t l1 = 0,l2 = 0;
      ierr = PetscStrlen(filename,&l1);CHKERRQ(ierr);
      ierr = PetscStrlen(imageext,&l2);CHKERRQ(ierr);
      ierr = PetscStrncpy(buf,filename,l1-l2+1);CHKERRQ(ierr);
      savename = buf;
    }
  }

  if (!savename) {ierr = PetscObjectGetName((PetscObject)draw,&savename);CHKERRQ(ierr);}
  ierr = PetscDrawImageCheckFormat(&imageext);CHKERRQ(ierr);
  if (movieext) {ierr = PetscDrawMovieCheckFormat(&movieext);CHKERRQ(ierr);}
  if (movieext) draw->savesinglefile = PETSC_FALSE; /* otherwise we cannot generage movies */

  if (draw->savesinglefile) {
    ierr = PetscInfo2(NULL,"Will save image to file %s%s\n",savename,imageext);CHKERRQ(ierr);
  } else {
    ierr = PetscInfo3(NULL,"Will save images to file %s/%s_%%d%s\n",savename,savename,imageext);CHKERRQ(ierr);
  }
  if (movieext) {
    ierr = PetscInfo2(NULL,"Will save movie to file %s%s\n",savename,movieext);CHKERRQ(ierr);
  }

  draw->savefilecount = 0;
  ierr = PetscFree(draw->savefilename);CHKERRQ(ierr);
  ierr = PetscFree(draw->saveimageext);CHKERRQ(ierr);
  ierr = PetscFree(draw->savemovieext);CHKERRQ(ierr);
  ierr = PetscStrallocpy(savename,&draw->savefilename);CHKERRQ(ierr);
  ierr = PetscStrallocpy(imageext,&draw->saveimageext);CHKERRQ(ierr);
  ierr = PetscStrallocpy(movieext,&draw->savemovieext);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*@C
   SNESMSRegister - register a multistage scheme

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

   Input Parameters:
+  name - identifier for method
.  nstages - number of stages
.  nregisters - number of registers used by low-storage implementation
.  gamma - coefficients, see Ketcheson's paper
.  delta - coefficients, see Ketcheson's paper
-  betasub - subdiagonal of Shu-Osher form

   Notes:
   The notation is described in Ketcheson (2010) Runge-Kutta methods with minimum storage implementations.

   Level: advanced

.keywords: SNES, register

.seealso: SNESMS
@*/
PetscErrorCode SNESMSRegister(const SNESMSType name,PetscInt nstages,PetscInt nregisters,PetscReal stability,const PetscReal gamma[],const PetscReal delta[],const PetscReal betasub[])
{
  PetscErrorCode ierr;
  SNESMSTableauLink link;
  SNESMSTableau t;

  PetscFunctionBegin;
  PetscValidCharPointer(name,1);
  if (nstages < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must have at least one stage");
  if (nregisters != 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only support for methods written in 3-register form");
  PetscValidPointer(gamma,4);
  PetscValidPointer(delta,5);
  PetscValidPointer(betasub,6);

  ierr = PetscMalloc(sizeof(*link),&link);CHKERRQ(ierr);
  ierr = PetscMemzero(link,sizeof(*link));CHKERRQ(ierr);
  t = &link->tab;
  ierr = PetscStrallocpy(name,&t->name);CHKERRQ(ierr);
  t->nstages    = nstages;
  t->nregisters = nregisters;
  t->stability  = stability;
  ierr = PetscMalloc3(nstages*nregisters,PetscReal,&t->gamma,nstages,PetscReal,&t->delta,nstages,PetscReal,&t->betasub);CHKERRQ(ierr);
  ierr = PetscMemcpy(t->gamma,gamma,nstages*nregisters*sizeof(PetscReal));CHKERRQ(ierr);
  ierr = PetscMemcpy(t->delta,delta,nstages*sizeof(PetscReal));CHKERRQ(ierr);
  ierr = PetscMemcpy(t->betasub,betasub,nstages*sizeof(PetscReal));CHKERRQ(ierr);

  link->next = SNESMSTableauList;
  SNESMSTableauList = link;
  PetscFunctionReturn(0);
}

// Logically collective
dErr dUnitsCreateUnit(dUnits un,const char *type,const char *longname,const char *shortname,dInt n,const dReal expon[],dUnit *newunit)
{
  dErr err;
  dUnit unit;

  dFunctionBegin;
  dValidHeader(un,dUNITS_CLASSID,1);
  if (n < 1 || n > dUNITS_MAX) dERROR(((dObject)un)->comm,PETSC_ERR_ARG_OUTOFRANGE,"The number of exponents %D must be positive, but no larger than %D",n,(dInt)dUNITS_MAX);
  dValidRealPointer(expon,5);
  dValidPointer(newunit,6);
  err = dUnitsGetEmptyUnit_Private(un,&unit);dCHK(err);
  err = PetscStrallocpy(type,&unit->quantity);dCHK(err);
  err = dUnitsAssignName(un,dUnitName,longname,n,expon,&unit->longname);dCHK(err);
  err = dUnitsAssignName(un,dUnitShortName,shortname,n,expon,&unit->shortname);dCHK(err);
  err = dUnitsAssignName(un,dUnitSIName,NULL,n,expon,&unit->siname);dCHK(err);
  unit->toSI = 1.0;
  unit->toCommon = 1.0;
  for (dInt i=0; i<n; i++) {
    dUnit base;
    err = dUnitsGetBase(un,i,&base);dCHK(err);
    unit->toCommon *= PetscPowScalar(dUnitDimensionalize(base,1.0),expon[i]);
    unit->toSI *= PetscPowScalar(dUnitDimensionalizeSI(base,1.0),expon[i]);
    unit->expon[i] = expon[i];
  }
  *newunit = unit;
  dFunctionReturn(0);
}

PetscErrorCode PetscViewerFileSetName_MPIIO(PetscViewer viewer,const char name[])
{
  PetscMPIInt        rank;
  PetscErrorCode     ierr;
  size_t             len;
  PetscViewer_Binary *vbinary = (PetscViewer_Binary*)viewer->data;
  char               *gz;
  PetscBool          found;
  PetscFileMode      type = vbinary->btype;

  PetscFunctionBegin;
  if (type == (PetscFileMode) -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call PetscViewerFileSetMode() before PetscViewerFileSetName()");
  ierr = PetscViewerFileClose_MPIIO(viewer);CHKERRQ(ierr);
  ierr = PetscOptionsGetBool(((PetscObject)viewer)->prefix,"-viewer_binary_skip_info",&vbinary->skipinfo,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetBool(((PetscObject)viewer)->prefix,"-viewer_binary_skip_options",&vbinary->skipoptions,NULL);CHKERRQ(ierr);

  ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)viewer),&rank);CHKERRQ(ierr);
  ierr = PetscStrallocpy(name,&vbinary->filename);CHKERRQ(ierr);

  vbinary->storecompressed = PETSC_FALSE;

  /* only first processor opens file if writeable */
  if (type == FILE_MODE_READ) {
    MPI_File_open(PetscObjectComm((PetscObject)viewer),vbinary->filename,MPI_MODE_RDONLY,MPI_INFO_NULL,&vbinary->mfdes);CHKERRQ(ierr);
  } else if (type == FILE_MODE_WRITE) {
    MPI_File_open(PetscObjectComm((PetscObject)viewer),vbinary->filename,MPI_MODE_WRONLY | MPI_MODE_CREATE,MPI_INFO_NULL,&vbinary->mfdes);CHKERRQ(ierr);
  }

  /*
      try to open info file: all processors open this file if read only

      Below is identical code to the code for Binary above, should be put in seperate routine
  */
  if (!vbinary->skipinfo && (!rank || type == FILE_MODE_READ)) {
    char infoname[PETSC_MAX_PATH_LEN],iname[PETSC_MAX_PATH_LEN];

    ierr = PetscStrcpy(infoname,name);CHKERRQ(ierr);
    /* remove .gz if it ends library name */
    ierr = PetscStrstr(infoname,".gz",&gz);CHKERRQ(ierr);
    if (gz) {
      ierr = PetscStrlen(gz,&len);CHKERRQ(ierr);
      if (len == 3) *gz = 0;
    }

    ierr = PetscStrcat(infoname,".info");CHKERRQ(ierr);
    ierr = PetscFixFilename(infoname,iname);CHKERRQ(ierr);
    if (type == FILE_MODE_READ) {
      ierr = PetscFileRetrieve(PetscObjectComm((PetscObject)viewer),iname,infoname,PETSC_MAX_PATH_LEN,&found);CHKERRQ(ierr);
      ierr = PetscOptionsInsertFile(PetscObjectComm((PetscObject)viewer),infoname,PETSC_FALSE);CHKERRQ(ierr);
    } else {
      vbinary->fdes_info = fopen(infoname,"w");
      if (!vbinary->fdes_info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN,"Cannot open .info file %s for writing",infoname);
    }
  }

#if defined(PETSC_USE_LOG)
  PetscLogObjectState((PetscObject)viewer,"File: %s",name);
#endif
  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);
}

PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc)
{
    PetscErrorCode ierr;
    PC_Cholesky    *dir;

    PetscFunctionBegin;
    ierr     = PetscNewLog(pc,&dir);
    CHKERRQ(ierr);
    pc->data = (void*)dir;
    ierr     = PCFactorInitialize(pc);
    CHKERRQ(ierr);

    ((PC_Factor*)dir)->factortype         = MAT_FACTOR_CHOLESKY;
    ((PC_Factor*)dir)->info.fill          = 5.0;

    dir->col = 0;
    dir->row = 0;

    ierr = PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)dir)->ordering);
    CHKERRQ(ierr);

    pc->ops->destroy           = PCDestroy_Cholesky;
    pc->ops->reset             = PCReset_Cholesky;
    pc->ops->apply             = PCApply_Cholesky;
    pc->ops->applytranspose    = PCApplyTranspose_Cholesky;
    pc->ops->setup             = PCSetUp_Cholesky;
    pc->ops->setfromoptions    = PCSetFromOptions_Cholesky;
    pc->ops->view              = PCView_Cholesky;
    pc->ops->applyrichardson   = 0;
    PetscFunctionReturn(0);
}

PetscErrorCode DMLabelDuplicate(DMLabel label, DMLabel *labelnew)
{
  PetscInt       v;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = DMLabelMakeAllValid_Private(label);CHKERRQ(ierr);
  ierr = PetscNew(labelnew);CHKERRQ(ierr);
  ierr = PetscStrallocpy(label->name, &(*labelnew)->name);CHKERRQ(ierr);

  (*labelnew)->refct        = 1;
  (*labelnew)->numStrata    = label->numStrata;
  (*labelnew)->defaultValue = label->defaultValue;
  if (label->numStrata) {
    ierr = PetscMalloc1(label->numStrata, &(*labelnew)->stratumValues);CHKERRQ(ierr);
    ierr = PetscMalloc1(label->numStrata, &(*labelnew)->stratumSizes);CHKERRQ(ierr);
    ierr = PetscMalloc1(label->numStrata, &(*labelnew)->ht);CHKERRQ(ierr);
    ierr = PetscMalloc1(label->numStrata, &(*labelnew)->points);CHKERRQ(ierr);
    ierr = PetscMalloc1(label->numStrata, &(*labelnew)->validIS);CHKERRQ(ierr);
    /* Could eliminate unused space here */
    for (v = 0; v < label->numStrata; ++v) {
      PetscHashICreate((*labelnew)->ht[v]);
      (*labelnew)->validIS[v]        = PETSC_TRUE;
      (*labelnew)->stratumValues[v]  = label->stratumValues[v];
      (*labelnew)->stratumSizes[v]   = label->stratumSizes[v];
      ierr = PetscObjectReference((PetscObject) (label->points[v]));CHKERRQ(ierr);
      (*labelnew)->points[v]         = label->points[v];
    }
  }
  (*labelnew)->pStart = -1;
  (*labelnew)->pEnd   = -1;
  (*labelnew)->bt     = NULL;
  PetscFunctionReturn(0);
}