PetscErrorCode VecCreate_Seq_Private(Vec v,const PetscScalar array[])
{
  Vec_Seq        *s;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = PetscNewLog(v,&s);CHKERRQ(ierr);
  ierr = PetscMemcpy(v->ops,&DvOps,sizeof(DvOps));CHKERRQ(ierr);

  v->data            = (void*)s;
  v->petscnative     = PETSC_TRUE;
  s->array           = (PetscScalar*)array;
  s->array_allocated = 0;

  ierr = PetscLayoutSetUp(v->map);CHKERRQ(ierr);
  ierr = PetscObjectChangeTypeName((PetscObject)v,VECSEQ);CHKERRQ(ierr);
#if defined(PETSC_HAVE_MATLAB_ENGINE)
  ierr = PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEnginePut_C",VecMatlabEnginePut_Default);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEngineGet_C",VecMatlabEngineGet_Default);CHKERRQ(ierr);
#endif
#if defined(PETSC_USE_MIXED_PRECISION)
  ((PetscObject)v)->precision = (PetscPrecision)sizeof(PetscReal);
#endif
  PetscFunctionReturn(0);
}

/*@C
   PFSetType - Builds PF for a particular function

   Collective on PF

   Input Parameter:
+  pf - the function context.
.  type - a known method
-  ctx - optional type dependent context

   Options Database Key:
.  -pf_type <type> - Sets PF type


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

  Level: intermediate

.keywords: PF, set, method, type

.seealso: PFSet(), PFRegisterDynamic(), PFCreate(), DMDACreatePF()

@*/
PetscErrorCode  PFSetType(PF pf,const PFType type,void *ctx)
{
  PetscErrorCode ierr,(*r)(PF,void*);
  PetscBool      match;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(pf,PF_CLASSID,1);
  PetscValidCharPointer(type,2);

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

  if (pf->ops->destroy) {ierr =  (*pf->ops->destroy)(pf);CHKERRQ(ierr);}
  pf->data        = 0;

  /* Determine the PFCreateXXX routine for a particular function */
  ierr =  PetscFListFind(PFList,((PetscObject)pf)->comm,type,PETSC_TRUE,(void (**)(void)) &r);CHKERRQ(ierr);
  if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested PF type %s",type);
  pf->ops->destroy             = 0;
  pf->ops->view                = 0;
  pf->ops->apply               = 0;
  pf->ops->applyvec            = 0;

  /* Call the PFCreateXXX routine for this particular function */
  ierr = (*r)(pf,ctx);CHKERRQ(ierr);

  ierr = PetscObjectChangeTypeName((PetscObject)pf,type);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode MatDestroy_DAAD(Mat A)
{
    Mat_DAAD       *a = (Mat_DAAD*)A->data;
    PetscErrorCode ierr;

    PetscFunctionBegin;
    ierr = DMDestroy(&a->da);
    CHKERRQ(ierr);
    ierr = VecDestroy(&a->localu);
    CHKERRQ(ierr);
    ierr = VecDestroy(&a->diagonal);
    CHKERRQ(ierr);
    ierr = PetscFree(A->data);
    CHKERRQ(ierr);
    ierr = PetscObjectChangeTypeName((PetscObject)A,0);
    CHKERRQ(ierr);
    ierr = PetscObjectComposeFunction((PetscObject)A,"MatMFFDSetBase_C","",PETSC_NULL);
    CHKERRQ(ierr);
    ierr = PetscObjectComposeFunction((PetscObject)A,"MatDAADSetDA_C","",PETSC_NULL);
    CHKERRQ(ierr);
    ierr = PetscObjectComposeFunction((PetscObject)A,"MatDAADSetSNES_C","",PETSC_NULL);
    CHKERRQ(ierr);
    ierr = PetscObjectComposeFunction((PetscObject)A,"MatDAADSetCtx_C","",PETSC_NULL);
    CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

/*@
      MatCreateTranspose - Creates a new matrix object that behaves like A'

   Collective on Mat

   Input Parameter:
.   A  - the (possibly rectangular) matrix

   Output Parameter:
.   N - the matrix that represents A'

   Level: intermediate

   Notes: The transpose A' is NOT actually formed! Rather the new matrix
          object performs the matrix-vector product by using the MatMultTranspose() on
          the original matrix

.seealso: MatCreateNormal(), MatMult(), MatMultTranspose(), MatCreate()

@*/
PetscErrorCode  MatCreateTranspose(Mat A,Mat *N)
{
  PetscErrorCode ierr;
  PetscInt       m,n;
  Mat_Transpose  *Na;

  PetscFunctionBegin;
  ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr);
  ierr = MatCreate(PetscObjectComm((PetscObject)A),N);CHKERRQ(ierr);
  ierr = MatSetSizes(*N,n,m,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
  ierr = PetscObjectChangeTypeName((PetscObject)*N,MATTRANSPOSEMAT);CHKERRQ(ierr);

  ierr       = PetscNewLog(*N,&Na);CHKERRQ(ierr);
  (*N)->data = (void*) Na;
  ierr       = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);
  Na->A      = A;

  (*N)->ops->destroy          = MatDestroy_Transpose;
  (*N)->ops->mult             = MatMult_Transpose;
  (*N)->ops->multadd          = MatMultAdd_Transpose;
  (*N)->ops->multtranspose    = MatMultTranspose_Transpose;
  (*N)->ops->multtransposeadd = MatMultTransposeAdd_Transpose;
  (*N)->assembled             = PETSC_TRUE;

  ierr = PetscLayoutSetBlockSize((*N)->rmap,A->cmap->bs);CHKERRQ(ierr);
  ierr = PetscLayoutSetBlockSize((*N)->cmap,A->rmap->bs);CHKERRQ(ierr);
  ierr = PetscLayoutSetUp((*N)->rmap);CHKERRQ(ierr);
  ierr = PetscLayoutSetUp((*N)->cmap);CHKERRQ(ierr);
  ierr = MatSetUp(*N);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "DMCreate_Composite"
PetscErrorCode  DMCreate_Composite(DM p)
{
  PetscErrorCode ierr;
  DM_Composite   *com;

  PetscFunctionBegin;
  ierr = PetscNewLog(p,DM_Composite,&com);CHKERRQ(ierr);
  p->data = com;
  ierr = PetscObjectChangeTypeName((PetscObject)p,"DMComposite");CHKERRQ(ierr);
  com->n            = 0;
  com->next         = PETSC_NULL;
  com->nDM          = 0;

  p->ops->createglobalvector              = DMCreateGlobalVector_Composite;
  p->ops->createlocalvector               = DMCreateLocalVector_Composite;
  p->ops->createlocaltoglobalmapping      = DMCreateLocalToGlobalMapping_Composite;
  p->ops->createlocaltoglobalmappingblock = 0;
  p->ops->createfieldis                   = DMCreateFieldIS_Composite;
  p->ops->createfielddecomposition        = DMCreateFieldDecomposition_Composite;
  p->ops->refine                          = DMRefine_Composite;
  p->ops->coarsen                         = DMCoarsen_Composite;
  p->ops->createinterpolation             = DMCreateInterpolation_Composite;
  p->ops->creatematrix                    = DMCreateMatrix_Composite;
  p->ops->getcoloring                     = DMCreateColoring_Composite;
  p->ops->globaltolocalbegin              = DMGlobalToLocalBegin_Composite;
  p->ops->globaltolocalend                = DMGlobalToLocalEnd_Composite;
  p->ops->destroy                         = DMDestroy_Composite;
  p->ops->view                            = DMView_Composite;
  p->ops->setup                           = DMSetUp_Composite;
  PetscFunctionReturn(0);
}

/*@C
   MatColoringSetType - Sets the type of coloring algorithm used

   Collective on MatColoring

   Input Parameter:
+  mc - the MatColoring context
-  type - the type of coloring

   Level: beginner

   Notes:  Possible types include the sequential types MATCOLORINGLF,
   MATCOLORINGSL, and MATCOLORINGID from the MINPACK package as well
   as a parallel MATCOLORINGMIS algorithm.

.keywords: Coloring, type

.seealso: MatColoringCreate(), MatColoringApply()
@*/
PetscErrorCode MatColoringSetType(MatColoring mc,MatColoringType type)
{
  PetscBool      match;
  PetscErrorCode ierr,(*r)(MatColoring);

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mc,MAT_COLORING_CLASSID,1);
  PetscValidCharPointer(type,2);
  ierr = PetscObjectTypeCompare((PetscObject)mc,type,&match);CHKERRQ(ierr);
  if (match) PetscFunctionReturn(0);
  ierr =  PetscFunctionListFind(MatColoringList,type,&r);CHKERRQ(ierr);
  if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested MatColoring type %s",type);
  if (mc->ops->destroy) {
    ierr             = (*(mc)->ops->destroy)(mc);CHKERRQ(ierr);
    mc->ops->destroy = NULL;
  }
  mc->ops->apply            = 0;
  mc->ops->view             = 0;
  mc->ops->setfromoptions   = 0;
  mc->ops->destroy          = 0;

  ierr = PetscObjectChangeTypeName((PetscObject)mc,type);CHKERRQ(ierr);
  ierr = (*r)(mc);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PETSC_EXTERN PetscErrorCode MatConvert_MPIAIJ_MPIAIJCRL(Mat A,MatType type,MatReuse reuse,Mat *newmat)
{
  PetscErrorCode ierr;
  Mat            B = *newmat;
  Mat_AIJCRL     *aijcrl;

  PetscFunctionBegin;
  if (reuse == MAT_INITIAL_MATRIX) {
    ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr);
  }

  ierr     = PetscNewLog(B,Mat_AIJCRL,&aijcrl);CHKERRQ(ierr);
  B->spptr = (void*) aijcrl;

  /* Set function pointers for methods that we inherit from AIJ but override. */
  B->ops->duplicate   = MatDuplicate_AIJCRL;
  B->ops->assemblyend = MatAssemblyEnd_MPIAIJCRL;
  B->ops->destroy     = MatDestroy_MPIAIJCRL;
  B->ops->mult        = MatMult_AIJCRL;

  /* If A has already been assembled, compute the permutation. */
  if (A->assembled) {
    ierr = MatMPIAIJCRL_create_aijcrl(B);CHKERRQ(ierr);
  }
  ierr    = PetscObjectChangeTypeName((PetscObject)B,MATMPIAIJCRL);CHKERRQ(ierr);
  *newmat = B;
  PetscFunctionReturn(0);
}

dErr VecCreateDohp(MPI_Comm comm,dInt bs,dInt n,dInt nc,dInt nghosts,const dInt ghosts[],Vec *v)
{
  Vec_MPI *vmpi;
  Vec      vc,vg;
  dScalar *a;
  dErr     err;

  dFunctionBegin;
  dValidPointer(v,7);
  *v = 0;
  err = VecCreateGhostBlock(comm,bs,nc*bs,PETSC_DECIDE,nghosts,ghosts,&vc);dCHK(err);
  err = VecGetArray(vc,&a);dCHK(err);
  err = VecCreateMPIWithArray(comm,n*bs,PETSC_DECIDE,a,&vg);dCHK(err);
  err = VecRestoreArray(vc,&a);dCHK(err);
  err = VecSetBlockSize(vg,bs);dCHK(err);
  vmpi = vg->data;
  if (vmpi->localrep) dERROR(PETSC_COMM_SELF,1,"Vector has localrep, expected no localrep");
  vmpi->localrep = vc;          /* subvert this field to mean closed rep */
  /* Since we subvect .localrep, VecDestroy_MPI will automatically destroy the closed form */
  vg->ops->duplicate = VecDuplicate_Dohp;
  //vg->ops->destroy   = VecDestroy_Dohp;
  /* It might be useful to set the (block) LocalToGlobal mapping here, but in the use case I have in mind, the user is
  * always working with the closed form anyway (in function evaluation).  The \e matrix does need a customized
  * LocalToGlobal mapping.
  */
  err = PetscObjectChangeTypeName((dObject)vg,VECDOHP);dCHK(err);
  *v = vg;
  dFunctionReturn(0);
}

PETSC_EXTERN PetscErrorCode DMCreate_Redundant(DM dm)
{
  PetscErrorCode ierr;
  DM_Redundant   *red;

  PetscFunctionBegin;
  ierr     = PetscNewLog(dm,&red);CHKERRQ(ierr);
  dm->data = red;

  ierr = PetscObjectChangeTypeName((PetscObject)dm,DMREDUNDANT);CHKERRQ(ierr);

  dm->ops->setup              = DMSetUp_Redundant;
  dm->ops->view               = DMView_Redundant;
  dm->ops->createglobalvector = DMCreateGlobalVector_Redundant;
  dm->ops->createlocalvector  = DMCreateLocalVector_Redundant;
  dm->ops->creatematrix       = DMCreateMatrix_Redundant;
  dm->ops->destroy            = DMDestroy_Redundant;
  dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Redundant;
  dm->ops->globaltolocalend   = DMGlobalToLocalEnd_Redundant;
  dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Redundant;
  dm->ops->localtoglobalend   = DMLocalToGlobalEnd_Redundant;
  dm->ops->refine             = DMRefine_Redundant;
  dm->ops->coarsen            = DMCoarsen_Redundant;
  dm->ops->createinterpolation= DMCreateInterpolation_Redundant;
  dm->ops->getcoloring        = DMCreateColoring_Redundant;

  ierr = PetscObjectComposeFunction((PetscObject)dm,"DMRedundantSetSize_C",DMRedundantSetSize_Redundant);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)dm,"DMRedundantGetSize_C",DMRedundantGetSize_Redundant);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*@C
   PetscViewerSetType - Builds PetscViewer for a particular implementation.

   Collective on PetscViewer

   Input Parameter:
+  viewer      - the PetscViewer context
-  type        - for example, PETSCVIEWERASCII

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

   Level: advanced

   Notes:
   See "include/petscviewer.h" for available methods (for instance,
   PETSCVIEWERSOCKET)

.seealso: PetscViewerCreate(), PetscViewerGetType(), PetscViewerType, PetscViewerSetFormat()
@*/
PetscErrorCode  PetscViewerSetType(PetscViewer viewer,PetscViewerType type)
{
  PetscErrorCode ierr,(*r)(PetscViewer);
  PetscBool      match;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,1);
  PetscValidCharPointer(type,2);
  ierr = PetscObjectTypeCompare((PetscObject)viewer,type,&match);CHKERRQ(ierr);
  if (match) PetscFunctionReturn(0);

  /* cleanup any old type that may be there */
  if (viewer->data) {
    ierr         = (*viewer->ops->destroy)(viewer);CHKERRQ(ierr);

    viewer->ops->destroy = NULL;
    viewer->data         = 0;
  }
  ierr = PetscMemzero(viewer->ops,sizeof(struct _PetscViewerOps));CHKERRQ(ierr);

  ierr =  PetscFunctionListFind(PetscViewerList,type,&r);CHKERRQ(ierr);
  if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown PetscViewer type given: %s",type);

  ierr = PetscObjectChangeTypeName((PetscObject)viewer,type);CHKERRQ(ierr);
  ierr = (*r)(viewer);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*@
  DMPlexClone - Creates a DMPlex object with the same mesh as the original.

  Collective on MPI_Comm

  Input Parameter:
. dm - The original DMPlex object

  Output Parameter:
. newdm  - The new DMPlex object

  Level: beginner

.keywords: DMPlex, create
@*/
PetscErrorCode DMPlexClone(DM dm, DM *newdm)
{
  DM_Plex        *mesh;
  Vec             coords;
  void           *ctx;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
  PetscValidPointer(newdm,2);
  ierr         = DMCreate(PetscObjectComm((PetscObject)dm), newdm);CHKERRQ(ierr);
  ierr         = PetscSFDestroy(&(*newdm)->sf);CHKERRQ(ierr);
  ierr         = PetscObjectReference((PetscObject) dm->sf);CHKERRQ(ierr);
  (*newdm)->sf = dm->sf;
  mesh         = (DM_Plex*) dm->data;
  mesh->refct++;
  (*newdm)->data = mesh;
  ierr           = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr);
  ierr           = DMInitialize_Plex(*newdm);CHKERRQ(ierr);
  ierr           = DMGetApplicationContext(dm, &ctx);CHKERRQ(ierr);
  ierr           = DMSetApplicationContext(*newdm, ctx);CHKERRQ(ierr);
  ierr           = DMGetCoordinatesLocal(dm, &coords);CHKERRQ(ierr);
  if (coords) {
    ierr         = DMSetCoordinatesLocal(*newdm, coords);CHKERRQ(ierr);
  } else {
    ierr         = DMGetCoordinates(dm, &coords);CHKERRQ(ierr);
    if (coords) {ierr = DMSetCoordinates(*newdm, coords);CHKERRQ(ierr);}
  }
  PetscFunctionReturn(0);
}

PetscErrorCode MatDestroy_MPIAdj(Mat mat)
{
  Mat_MPIAdj     *a = (Mat_MPIAdj*)mat->data;
  PetscErrorCode ierr;

  PetscFunctionBegin;
#if defined(PETSC_USE_LOG)
  PetscLogObjectState((PetscObject)mat,"Rows=%D, Cols=%D, NZ=%D",mat->rmap->n,mat->cmap->n,a->nz);
#endif
  ierr = PetscFree(a->diag);CHKERRQ(ierr);
  if (a->freeaij) {
    if (a->freeaijwithfree) {
      if (a->i) free(a->i);
      if (a->j) free(a->j);
    } else {
      ierr = PetscFree(a->i);CHKERRQ(ierr);
      ierr = PetscFree(a->j);CHKERRQ(ierr);
      ierr = PetscFree(a->values);CHKERRQ(ierr);
    }
  }
  ierr = PetscFree(a->rowvalues);CHKERRQ(ierr);
  ierr = PetscFree(mat->data);CHKERRQ(ierr);
  ierr = PetscObjectChangeTypeName((PetscObject)mat,0);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMPIAdjSetPreallocation_C",NULL);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMPIAdjCreateNonemptySubcommMat_C",NULL);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*@C
    MatMFFDSetType - Sets the method that is used to compute the
    differencing parameter for finite differene matrix-free formulations.

    Input Parameters:
+   mat - the "matrix-free" matrix created via MatCreateSNESMF(), or MatCreateMFFD()
          or MatSetType(mat,MATMFFD);
-   ftype - the type requested, either MATMFFD_WP or MATMFFD_DS

    Level: advanced

    Notes:
    For example, such routines can compute h for use in
    Jacobian-vector products of the form

                        F(x+ha) - F(x)
          F'(u)a  ~=  ----------------
                              h

.seealso: MatCreateSNESMF(), MatMFFDRegister(), MatMFFDSetFunction()
@*/
PetscErrorCode  MatMFFDSetType(Mat mat,MatMFFDType ftype)
{
  PetscErrorCode ierr,(*r)(MatMFFD);
  MatMFFD        ctx = (MatMFFD)mat->data;
  PetscBool      match;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
  PetscValidCharPointer(ftype,2);

  ierr = PetscObjectTypeCompare((PetscObject)mat,MATMFFD,&match);CHKERRQ(ierr);
  if (!match) PetscFunctionReturn(0);

  /* already set, so just return */
  ierr = PetscObjectTypeCompare((PetscObject)ctx,ftype,&match);CHKERRQ(ierr);
  if (match) PetscFunctionReturn(0);

  /* destroy the old one if it exists */
  if (ctx->ops->destroy) {
    ierr = (*ctx->ops->destroy)(ctx);CHKERRQ(ierr);
  }

  ierr =  PetscFunctionListFind(MatMFFDList,ftype,&r);CHKERRQ(ierr);
  if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown MatMFFD type %s given",ftype);
  ierr = (*r)(ctx);CHKERRQ(ierr);
  ierr = PetscObjectChangeTypeName((PetscObject)ctx,ftype);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PETSC_EXTERN PetscErrorCode DMCreate_Shell(DM dm)
{
    PetscErrorCode ierr;
    DM_Shell       *shell;

    PetscFunctionBegin;
    ierr     = PetscNewLog(dm,&shell);
    CHKERRQ(ierr);
    dm->data = shell;

    ierr = PetscObjectChangeTypeName((PetscObject)dm,DMSHELL);
    CHKERRQ(ierr);

    dm->ops->destroy            = DMDestroy_Shell;
    dm->ops->createglobalvector = DMCreateGlobalVector_Shell;
    dm->ops->createlocalvector  = DMCreateLocalVector_Shell;
    dm->ops->creatematrix       = DMCreateMatrix_Shell;
    dm->ops->view               = DMView_Shell;
    dm->ops->load               = DMLoad_Shell;
    dm->ops->globaltolocalbegin = DMGlobalToLocalBeginDefaultShell;
    dm->ops->globaltolocalend   = DMGlobalToLocalEndDefaultShell;
    dm->ops->localtoglobalbegin = DMLocalToGlobalBeginDefaultShell;
    dm->ops->localtoglobalend   = DMLocalToGlobalEndDefaultShell;
    dm->ops->localtolocalbegin  = DMLocalToLocalBeginDefaultShell;
    dm->ops->localtolocalend    = DMLocalToLocalEndDefaultShell;
    dm->ops->createsubdm        = DMCreateSubDM_Shell;
    PetscFunctionReturn(0);
}

/*@C
   PCSetType - Builds PC for a particular preconditioner type

   Collective on PC

   Input Parameter:
+  pc - the preconditioner context.
-  type - a known method

   Options Database Key:
.  -pc_type <type> - Sets PC type

   Use -help for a list of available methods (for instance,
   jacobi or bjacobi)

  Notes:
  See "petsc/include/petscpc.h" for available methods (for instance,
  PCJACOBI, PCILU, or PCBJACOBI).

  Normally, it is best to use the KSPSetFromOptions() command and
  then set the PC type from the options database rather than by using
  this routine.  Using the options database provides the user with
  maximum flexibility in evaluating the many different preconditioners.
  The PCSetType() routine is provided for those situations where it
  is necessary to set the preconditioner independently of the command
  line or options database.  This might be the case, for example, when
  the choice of preconditioner changes during the execution of the
  program, and the user's application is taking responsibility for
  choosing the appropriate preconditioner.  In other words, this
  routine is not for beginners.

  Level: intermediate

  Developer Note: PCRegister() is used to add preconditioner types to PCList from which they
  are accessed by PCSetType().

.keywords: PC, set, method, type

.seealso: KSPSetType(), PCType, PCRegister(), PCCreate(), KSPGetPC()

@*/
PetscErrorCode  PCSetType(PC pc,PCType type)
{
  PetscErrorCode ierr,(*r)(PC);
  PetscBool      match;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc,PC_CLASSID,1);
  PetscValidCharPointer(type,2);

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

  ierr =  PetscFunctionListFind(PCList,type,&r);CHKERRQ(ierr);
  if (!r) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested PC type %s",type);
  /* Destroy the previous private PC context */
  if (pc->ops->destroy) {
    ierr             =  (*pc->ops->destroy)(pc);CHKERRQ(ierr);
    pc->ops->destroy = NULL;
    pc->data         = 0;
  }
  ierr = PetscFunctionListDestroy(&((PetscObject)pc)->qlist);CHKERRQ(ierr);
  /* Reinitialize function pointers in PCOps structure */
  ierr = PetscMemzero(pc->ops,sizeof(struct _PCOps));CHKERRQ(ierr);
  /* XXX Is this OK?? */
  pc->modifysubmatrices  = 0;
  pc->modifysubmatricesP = 0;
  /* Call the PCCreate_XXX routine for this particular preconditioner */
  pc->setupcalled = 0;

  ierr = PetscObjectChangeTypeName((PetscObject)pc,type);CHKERRQ(ierr);
  ierr = (*r)(pc);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*@
      MatCreateHermitianTranspose - Creates a new matrix object that behaves like A'*

   Collective on Mat

   Input Parameter:
.   A  - the (possibly rectangular) matrix

   Output Parameter:
.   N - the matrix that represents A'*

   Level: intermediate

   Notes:
    The hermitian transpose A' is NOT actually formed! Rather the new matrix
          object performs the matrix-vector product by using the MatMultHermitianTranspose() on
          the original matrix

.seealso: MatCreateNormal(), MatMult(), MatMultHermitianTranspose(), MatCreate()

@*/
PetscErrorCode  MatCreateHermitianTranspose(Mat A,Mat *N)
{
  PetscErrorCode ierr;
  PetscInt       m,n;
  Mat_HT         *Na;

  PetscFunctionBegin;
  ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr);
  ierr = MatCreate(PetscObjectComm((PetscObject)A),N);CHKERRQ(ierr);
  ierr = MatSetSizes(*N,n,m,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
  ierr = PetscLayoutSetUp((*N)->rmap);CHKERRQ(ierr);
  ierr = PetscLayoutSetUp((*N)->cmap);CHKERRQ(ierr);
  ierr = PetscObjectChangeTypeName((PetscObject)*N,MATTRANSPOSEMAT);CHKERRQ(ierr);

  ierr       = PetscNewLog(*N,&Na);CHKERRQ(ierr);
  (*N)->data = (void*) Na;
  ierr       = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);
  Na->A      = A;

  (*N)->ops->destroy                    = MatDestroy_HT;
  (*N)->ops->mult                       = MatMult_HT;
  (*N)->ops->multadd                    = MatMultAdd_HT;
  (*N)->ops->multhermitiantranspose     = MatMultHermitianTranspose_HT;
  (*N)->ops->multhermitiantransposeadd  = MatMultHermitianTransposeAdd_HT;
  (*N)->ops->duplicate                  = MatDuplicate_HT;
  (*N)->ops->getvecs                    = MatCreateVecs_HT;
  (*N)->ops->axpy                       = MatAXPY_HT;
  (*N)->assembled                       = PETSC_TRUE;

  ierr = PetscObjectComposeFunction((PetscObject)(*N),"MatHermitianTransposeGetMat_C",MatHermitianTransposeGetMat_HT);CHKERRQ(ierr);
  ierr = MatSetBlockSizes(*N,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr);
  ierr = MatSetUp(*N);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PETSC_EXTERN PetscErrorCode MatConvert_MPISBAIJ_MPISBSTRM(Mat A,MatType type,MatReuse reuse,Mat *newmat)
{
  PetscErrorCode ierr;
  Mat            B = *newmat;
  Mat_SeqSBSTRM  *sbstrm;

  PetscFunctionBegin;
  if (reuse == MAT_INITIAL_MATRIX) {
    ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr);
  }
  /* printf(" --- in MatConvert_MPISBAIJ_MPISBSTRM  -- 1 \n"); */

  ierr     = PetscNewLog(B,   Mat_SeqSBSTRM,&sbstrm);CHKERRQ(ierr);
  B->spptr = (void*)sbstrm;

  /* Set function pointers for methods that we inherit from AIJ but override.
     B->ops->duplicate   = MatDuplicate_SBSTRM;
     B->ops->mult        = MatMult_SBSTRM;
     B->ops->destroy     = MatDestroy_MPISBSTRM;
  */
  B->ops->assemblyend = MatAssemblyEnd_MPISBSTRM;

  /* If A has already been assembled, compute the permutation. */
  if (A->assembled) {
    ierr = MPISBSTRM_create_sbstrm(B);CHKERRQ(ierr);
  }

  ierr = PetscObjectChangeTypeName((PetscObject) B, MATMPISBSTRM);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)B,"MatMPISBAIJSetPreallocation_C",MatMPISBAIJSetPreallocation_MPISBSTRM);CHKERRQ(ierr);
  *newmat = B;
  PetscFunctionReturn(0);
}

/*@C
   EPSSetType - Selects the particular solver to be used in the EPS object.

   Logically Collective on EPS

   Input Parameters:
+  eps  - the eigensolver context
-  type - a known method

   Options Database Key:
.  -eps_type <method> - Sets the method; use -help for a list
    of available methods

   Notes:
   See "slepc/include/slepceps.h" for available methods. The default
   is EPSKRYLOVSCHUR.

   Normally, it is best to use the EPSSetFromOptions() command and
   then set the EPS type from the options database rather than by using
   this routine.  Using the options database provides the user with
   maximum flexibility in evaluating the different available methods.
   The EPSSetType() routine is provided for those situations where it
   is necessary to set the iterative solver independently of the command
   line or options database.

   Level: intermediate

.seealso: STSetType(), EPSType
@*/
PetscErrorCode EPSSetType(EPS eps,EPSType type)
{
    PetscErrorCode ierr,(*r)(EPS);
    PetscBool      match;

    PetscFunctionBegin;
    PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
    PetscValidCharPointer(type,2);

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

    ierr = PetscFunctionListFind(EPSList,type,&r);
    CHKERRQ(ierr);
    if (!r) SETERRQ1(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown EPS type given: %s",type);

    if (eps->ops->destroy) {
        ierr = (*eps->ops->destroy)(eps);
        CHKERRQ(ierr);
    }
    ierr = PetscMemzero(eps->ops,sizeof(struct _EPSOps));
    CHKERRQ(ierr);

    eps->state = EPS_STATE_INITIAL;
    ierr = PetscObjectChangeTypeName((PetscObject)eps,type);
    CHKERRQ(ierr);
    ierr = (*r)(eps);
    CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

static PetscErrorCode PCDestroy_PARMS(PC pc)
{
  PC_PARMS       *parms = (PC_PARMS*)pc->data;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  if (parms->map) parms_MapFree(&parms->map);
  if (parms->A) parms_MatFree(&parms->A);
  if (parms->pc) parms_PCFree(&parms->pc);
  if (parms->lvec0) {
    ierr = PetscFree(parms->lvec0);CHKERRQ(ierr);
  }
  if (parms->lvec1) {
    ierr = PetscFree(parms->lvec1);CHKERRQ(ierr);
  }
  ierr = PetscFree(pc->data);CHKERRQ(ierr);

  ierr = PetscObjectChangeTypeName((PetscObject)pc,0);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)pc,"PCPARMSSetGlobal_C",NULL);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)pc,"PCPARMSSetLocal_C",NULL);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)pc,"PCPARMSSetSolveTolerances_C",NULL);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)pc,"PCPARMSSetSolveRestart_C",NULL);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)pc,"PCPARMSSetNonsymPerm_C",NULL);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)pc,"PCPARMSSetFill_C",NULL);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PETSC_EXTERN PetscErrorCode MatCreate_Preallocator(Mat A)
{
  Mat_Preallocator *p;
  PetscErrorCode    ierr;

  PetscFunctionBegin;
  ierr = PetscNewLog(A, &p);CHKERRQ(ierr);
  A->data = (void *) p;

  p->ht  = NULL;
  p->dnz = NULL;
  p->onz = NULL;

  /* matrix ops */
  ierr = PetscMemzero(A->ops, sizeof(struct _MatOps));CHKERRQ(ierr);
  A->ops->destroy                 = MatDestroy_Preallocator;
  A->ops->setup                   = MatSetUp_Preallocator;
  A->ops->setvalues               = MatSetValues_Preallocator;
  A->ops->assemblybegin           = MatAssemblyBegin_Preallocator;
  A->ops->assemblyend             = MatAssemblyEnd_Preallocator;
  A->ops->view                    = MatView_Preallocator;
  A->ops->setoption               = MatSetOption_Preallocator;

  /* special MATPREALLOCATOR functions */
  ierr = PetscObjectComposeFunction((PetscObject) A, "MatPreallocatorPreallocate_C", MatPreallocatorPreallocate_Preallocator);CHKERRQ(ierr);
  ierr = PetscObjectChangeTypeName((PetscObject) A, MATPREALLOCATOR);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}