EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "PCCreate_CP"
PetscErrorCode  PCCreate_CP(PC pc)
{
  PC_CP          *cp;
  PetscErrorCode ierr;

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

  pc->ops->apply               = PCApply_CP;
  pc->ops->applytranspose      = PCApply_CP;
  pc->ops->setup               = PCSetUp_CP;
  pc->ops->reset               = PCReset_CP;
  pc->ops->destroy             = PCDestroy_CP;
  pc->ops->setfromoptions      = PCSetFromOptions_CP;
  pc->ops->view                = 0;
  pc->ops->applyrichardson     = 0;
  PetscFunctionReturn(0);
}

PETSC_EXTERN PetscErrorCode MatPartitioningCreate_Hierarchical(MatPartitioning part)
{
  PetscErrorCode                ierr;
  MatPartitioning_Hierarchical *hpart;

  PetscFunctionBegin;
  ierr       = PetscNewLog(part,&hpart);CHKERRQ(ierr);
  part->data = (void*)hpart;

  hpart->fineparttype       = 0; /* fine level partitioner */
  hpart->coarseparttype     = 0; /* coarse level partitioner */
  hpart->Nfineparts         = 1; /* we do not further partition coarse partition any more by default */
  hpart->Ncoarseparts       = 0; /* number of coarse parts (first level) */
  hpart->coarseparts        = 0;
  hpart->fineparts          = 0;

  part->ops->apply          = MatPartitioningApply_Hierarchical;
  part->ops->view           = MatPartitioningView_Hierarchical;
  part->ops->destroy        = MatPartitioningDestroy_Hierarchical;
  part->ops->setfromoptions = MatPartitioningSetFromOptions_Hierarchical;
  PetscFunctionReturn(0);
}

PETSC_EXTERN PetscErrorCode PCCreate_KSP(PC pc)
{
  PetscErrorCode ierr;
  PC_KSP         *jac;

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

  ierr = PetscMemzero(pc->ops,sizeof(struct _PCOps));CHKERRQ(ierr);
  pc->ops->apply           = PCApply_KSP;
  pc->ops->applytranspose  = PCApplyTranspose_KSP;
  pc->ops->setup           = PCSetUp_KSP;
  pc->ops->reset           = PCReset_KSP;
  pc->ops->destroy         = PCDestroy_KSP;
  pc->ops->setfromoptions  = PCSetFromOptions_KSP;
  pc->ops->view            = PCView_KSP;

  ierr = PetscObjectComposeFunction((PetscObject)pc,"PCKSPGetKSP_C",PCKSPGetKSP_KSP);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)pc,"PCKSPSetKSP_C",PCKSPSetKSP_KSP);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PETSC_EXTERN PetscErrorCode NEPCreate_Interpol(NEP nep)
{
  PetscErrorCode ierr;
  NEP_INTERPOL   *ctx;

  PetscFunctionBegin;
  ierr = PetscNewLog(nep,&ctx);CHKERRQ(ierr);
  ctx->deg  = 5;
  nep->data = (void*)ctx;

  nep->ops->solve          = NEPSolve_Interpol;
  nep->ops->setup          = NEPSetUp_Interpol;
  nep->ops->setfromoptions = NEPSetFromOptions_Interpol;
  nep->ops->reset          = NEPReset_Interpol;
  nep->ops->destroy        = NEPDestroy_Interpol;
  nep->ops->view           = NEPView_Interpol;
  ierr = PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolSetDegree_C",NEPInterpolSetDegree_Interpol);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolGetDegree_C",NEPInterpolGetDegree_Interpol);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolSetPEP_C",NEPInterpolSetPEP_Interpol);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolGetPEP_C",NEPInterpolGetPEP_Interpol);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*@
      MatCreateFFT - Creates a matrix object that provides FFT via an external package

   Collective on MPI_Comm

   Input Parameter:
+   comm - MPI communicator
.   ndim - the ndim-dimensional transform
.   dim - array of size ndim, dim[i] contains the vector length in the i-dimension
-   type - package type, e.g., FFTW or FFTCU

   Output Parameter:
.   A  - the matrix

  Options Database Keys:
+ -mat_fft_type - set FFT type

   Level: intermediate

@*/
PetscErrorCode MatCreateFFT(MPI_Comm comm,PetscInt ndim,const PetscInt dim[],MatType mattype,Mat* A)
{
  PetscErrorCode ierr;
  PetscMPIInt    size;
  Mat            FFT;
  PetscInt       N,i;
  Mat_FFT        *fft;

  PetscFunctionBegin;
  if (ndim < 1) SETERRQ1(comm,PETSC_ERR_USER,"ndim %d must be > 0",ndim);
  ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);

  ierr = MatCreate(comm,&FFT);CHKERRQ(ierr);
  ierr = PetscNewLog(FFT,Mat_FFT,&fft);CHKERRQ(ierr);
  FFT->data = (void*)fft;
  N = 1;
  for (i=0; i<ndim; i++){
    if (dim[i] < 1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"dim[%d]=%d must be > 0",i,dim[i]);
    N *= dim[i];
  }

  ierr = PetscMalloc(ndim*sizeof(PetscInt),&fft->dim);CHKERRQ(ierr);
  ierr = PetscMemcpy(fft->dim,dim,ndim*sizeof(PetscInt));CHKERRQ(ierr);

  fft->ndim = ndim;
  fft->n    = PETSC_DECIDE;
  fft->N    = N;
  fft->data = PETSC_NULL;

  ierr = MatSetType(FFT,mattype);CHKERRQ(ierr);
  FFT->ops->destroy = MatDestroy_FFT;

  /* get runtime options */
  ierr = PetscOptionsBegin(((PetscObject)FFT)->comm,((PetscObject)FFT)->prefix,"FFT Options","Mat");CHKERRQ(ierr);
  PetscOptionsEnd();

  *A = FFT;
  PetscFunctionReturn(0);
}

EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "TSCreate_RK"
PetscErrorCode  TSCreate_RK(TS ts)
{
  TS_RK          *rk;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ts->ops->setup           = TSSetUp_RK;
  ts->ops->solve           = TSSolve_RK;
  ts->ops->destroy         = TSDestroy_RK;
  ts->ops->setfromoptions  = TSSetFromOptions_RK;
  ts->ops->view            = TSView_RK;

  ierr = PetscNewLog(ts,TS_RK,&rk);CHKERRQ(ierr);
  ts->data = (void*)rk;

  ierr = PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSRKSetTolerance_C","TSRKSetTolerance_RK",TSRKSetTolerance_RK);CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

/*MC
   SNESLINESEARCHBT - Backtracking line search.

   This line search finds the minimum of a polynomial fitting of the L2 norm of the
   function. If this fit does not satisfy the conditions for progress, the interval shrinks
   and the fit is reattempted at most max_it times or until lambda is below minlambda.

   Options Database Keys:
+  -snes_linesearch_alpha<1e-4> - slope descent parameter
.  -snes_linesearch_damping<1.0> - initial step length
.  -snes_linesearch_max_it<40> - maximum number of shrinking step
.  -snes_linesearch_minlambda<1e-12> - minimum step length allowed
-  -snes_linesearch_order<cubic,quadratic> - order of the approximation

   Level: advanced

   Notes:
   This line search is taken from "Numerical Methods for Unconstrained
   Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325.

.keywords: SNES, SNESLineSearch, damping

.seealso: SNESLineSearchCreate(), SNESLineSearchSetType()
M*/
PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_BT(SNESLineSearch linesearch)
{

  SNESLineSearch_BT *bt;
  PetscErrorCode    ierr;

  PetscFunctionBegin;
  linesearch->ops->apply          = SNESLineSearchApply_BT;
  linesearch->ops->destroy        = SNESLineSearchDestroy_BT;
  linesearch->ops->setfromoptions = SNESLineSearchSetFromOptions_BT;
  linesearch->ops->reset          = NULL;
  linesearch->ops->view           = SNESLineSearchView_BT;
  linesearch->ops->setup          = NULL;

  ierr = PetscNewLog(linesearch, SNESLineSearch_BT, &bt);CHKERRQ(ierr);

  linesearch->data    = (void*)bt;
  linesearch->max_its = 40;
  linesearch->order   = SNES_LINESEARCH_ORDER_CUBIC;
  bt->alpha           = 1e-4;
  PetscFunctionReturn(0);
}

EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "DMCreate_Cartesian"
PetscErrorCode DMCreate_Cartesian(DM dm)
{
    DM_Cartesian  *mesh;
    PetscErrorCode ierr;

    PetscFunctionBegin;
    PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
    ierr = PetscNewLog(dm, DM_Cartesian, &mesh);
    CHKERRQ(ierr);
    dm->data = mesh;

    new(&mesh->m) ALE::Obj<ALE::CartesianMesh>(PETSC_NULL);

    ierr = PetscStrallocpy(VECSTANDARD, &dm->vectype);
    CHKERRQ(ierr);
    dm->ops->globaltolocalbegin = 0;
    dm->ops->globaltolocalend   = 0;
    dm->ops->localtoglobalbegin = 0;
    dm->ops->localtoglobalend   = 0;
    dm->ops->createglobalvector = 0; /* DMCreateGlobalVector_Cartesian; */
    dm->ops->createlocalvector  = 0; /* DMCreateLocalVector_Cartesian; */
    dm->ops->createinterpolation   = DMCreateInterpolation_Cartesian;
    dm->ops->getcoloring        = 0;
    dm->ops->creatematrix          = 0; /* DMCreateMatrix_Cartesian; */
    dm->ops->refine             = DMRefine_Cartesian;
    dm->ops->coarsen            = DMCoarsen_Cartesian;
    dm->ops->refinehierarchy    = 0;
    dm->ops->coarsenhierarchy   = 0;
    dm->ops->getinjection       = 0;
    dm->ops->getaggregates      = 0;
    dm->ops->destroy            = DMDestroy_Cartesian;
    dm->ops->view               = DMView_Cartesian;
    dm->ops->setfromoptions     = DMSetFromOptions_Cartesian;
    dm->ops->setup              = 0;
    PetscFunctionReturn(0);
}

PetscErrorCode SNESDiffParameterCreate_More(SNES snes,Vec x,void **outneP)
{
  DIFFPAR_MORE   *neP;
  Vec            w;
  PetscRandom    rctx;  /* random number generator context */
  PetscErrorCode ierr;
  PetscBool      flg;
  char           noise_file[PETSC_MAX_PATH_LEN];

  PetscFunctionBegin;
  ierr = PetscNewLog(snes,&neP);CHKERRQ(ierr);

  neP->function_count = 0;
  neP->fnoise_min     = 1.0e-20;
  neP->hopt_min       = 1.0e-8;
  neP->h_first_try    = 1.0e-3;
  neP->fnoise_resets  = 0;
  neP->hopt_resets    = 0;

  /* Create work vectors */
  ierr = VecDuplicateVecs(x,3,&neP->workv);CHKERRQ(ierr);
  w    = neP->workv[0];

  /* Set components of vector w to random numbers */
  ierr = PetscRandomCreate(PetscObjectComm((PetscObject)snes),&rctx);CHKERRQ(ierr);
  ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
  ierr = VecSetRandom(w,rctx);CHKERRQ(ierr);
  ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr);

  /* Open output file */
  ierr = PetscOptionsGetString(((PetscObject)snes)->prefix,"-snes_mf_noise_file",noise_file,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
  if (flg) neP->fp = fopen(noise_file,"w");
  else     neP->fp = fopen("noise.out","w");
  if (!neP->fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN,"Cannot open file");
  ierr = PetscInfo(snes,"Creating Jorge's differencing parameter context\n");CHKERRQ(ierr);

  *outneP = neP;
  PetscFunctionReturn(0);
}

/*MC
      SNESNEWTONLS - Newton based nonlinear solver that uses a line search

   Options Database:
+   -snes_linesearch_type <bt> - bt,basic.  Select line search type
.   -snes_linesearch_order <3> - 2, 3. Selects the order of the line search for bt
.   -snes_linesearch_norms <true> - Turns on/off computation of the norms for basic linesearch
.   -snes_linesearch_alpha <alpha> - Sets alpha used in determining if reduction in function norm is sufficient
.   -snes_linesearch_maxstep <maxstep> - Sets the maximum stepsize the line search will use (if the 2-norm(y) > maxstep then scale y to be y = (maxstep/2-norm(y)) *y)
.   -snes_linesearch_minlambda <minlambda>  - Sets the minimum lambda the line search will tolerate
.   -snes_linesearch_monitor - print information about progress of line searches
-   -snes_linesearch_damping - damping factor used for basic line search

    Notes: This is the default nonlinear solver in SNES

   Level: beginner

.seealso:  SNESCreate(), SNES, SNESSetType(), SNESNEWTONTR, SNESQN, SNESLineSearchSetType(), SNESLineSearchSetOrder()
           SNESLineSearchSetPostCheck(), SNESLineSearchSetPreCheck() SNESLineSearchSetComputeNorms()

M*/
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "SNESCreate_NEWTONLS"
PetscErrorCode  SNESCreate_NEWTONLS(SNES snes)
{
  PetscErrorCode ierr;
  SNES_NEWTONLS  *neP;

  PetscFunctionBegin;
  snes->ops->setup           = SNESSetUp_NEWTONLS;
  snes->ops->solve           = SNESSolve_NEWTONLS;
  snes->ops->destroy         = SNESDestroy_NEWTONLS;
  snes->ops->setfromoptions  = SNESSetFromOptions_NEWTONLS;
  snes->ops->view            = SNESView_NEWTONLS;
  snes->ops->reset           = SNESReset_NEWTONLS;

  snes->usesksp                      = PETSC_TRUE;
  snes->usespc                       = PETSC_FALSE;
  ierr                               = PetscNewLog(snes,SNES_NEWTONLS,&neP);CHKERRQ(ierr);
  snes->data                         = (void*)neP;
  PetscFunctionReturn(0);
}

PETSC_EXTERN PetscErrorCode KSPCreate_PIPEFGMRES(KSP ksp)
{
  KSP_PIPEFGMRES *pipefgmres;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = PetscNewLog(ksp,&pipefgmres);CHKERRQ(ierr);

  ksp->data                              = (void*)pipefgmres;
  ksp->ops->buildsolution                = KSPBuildSolution_PIPEFGMRES;
  ksp->ops->setup                        = KSPSetUp_PIPEFGMRES;
  ksp->ops->solve                        = KSPSolve_PIPEFGMRES;
  ksp->ops->reset                        = KSPReset_PIPEFGMRES;
  ksp->ops->destroy                      = KSPDestroy_PIPEFGMRES;
  ksp->ops->view                         = KSPView_PIPEFGMRES;
  ksp->ops->setfromoptions               = KSPSetFromOptions_PIPEFGMRES;
  ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_GMRES;
  ksp->ops->computeeigenvalues           = KSPComputeEigenvalues_GMRES;

  ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,3);CHKERRQ(ierr);

  ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetPreAllocateVectors_C",KSPGMRESSetPreAllocateVectors_GMRES);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetRestart_C",KSPGMRESSetRestart_GMRES);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetRestart_C",KSPGMRESGetRestart_GMRES);CHKERRQ(ierr);

  pipefgmres->nextra_vecs    = 1;
  pipefgmres->haptol         = 1.0e-30;
  pipefgmres->q_preallocate  = 0;
  pipefgmres->delta_allocate = PIPEFGMRES_DELTA_DIRECTIONS;
  pipefgmres->orthog         = 0;
  pipefgmres->nrs            = 0;
  pipefgmres->sol_temp       = 0;
  pipefgmres->max_k          = PIPEFGMRES_DEFAULT_MAXK;
  pipefgmres->Rsvd           = 0;
  pipefgmres->orthogwork     = 0;
  pipefgmres->cgstype        = KSP_GMRES_CGS_REFINE_NEVER;
  pipefgmres->shift          = 1.0;
  PetscFunctionReturn(0);
}

/*@C
   SNESMonitorSetRatio - Sets SNES to use a monitor that prints the
   ratio of the function norm at each iteration.

   Collective on SNES

   Input Parameters:
+   snes - the SNES context
-   viewer - ASCII viewer to print output

   Level: intermediate

.keywords: SNES, nonlinear, monitor, norm

.seealso: SNESMonitorSet(), SNESMonitorSolution(), SNESMonitorDefault()
@*/
PetscErrorCode  SNESMonitorSetRatio(SNES snes,PetscViewer viewer)
{
  PetscErrorCode          ierr;
  SNESMonitorRatioContext *ctx;
  PetscReal               *history;

  PetscFunctionBegin;
  if (!viewer) {
    ierr = PetscViewerASCIIOpen(PetscObjectComm((PetscObject)snes),"stdout",&viewer);CHKERRQ(ierr);
    ierr = PetscObjectReference((PetscObject)viewer);CHKERRQ(ierr);
  }
  ierr = PetscNewLog(snes,SNESMonitorRatioContext,&ctx);CHKERRQ(ierr);
  ierr = SNESGetConvergenceHistory(snes,&history,NULL,NULL);CHKERRQ(ierr);
  if (!history) {
    ierr = PetscMalloc(100*sizeof(PetscReal),&ctx->history);CHKERRQ(ierr);
    ierr = SNESSetConvergenceHistory(snes,ctx->history,0,100,PETSC_TRUE);CHKERRQ(ierr);
  }
  ctx->viewer = viewer;

  ierr = SNESMonitorSet(snes,SNESMonitorRatio,ctx,SNESMonitorRatioDestroy);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat A,MatFactorType ftype,Mat *F)
{
  Mat            B;
  Mat_CHOLMOD    *chol;
  PetscErrorCode ierr;
  PetscInt       m=A->rmap->n,n=A->cmap->n,bs;

  PetscFunctionBegin;
  if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"CHOLMOD cannot do %s factorization with SBAIJ, only %s",
                                             MatFactorTypes[ftype],MatFactorTypes[MAT_FACTOR_CHOLESKY]);
  ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr);
  if (bs != 1) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"CHOLMOD only supports block size=1, given %D",bs);
  /* Create the factorization matrix F */
  ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr);
  ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr);
  ierr = PetscStrallocpy("cholmod",&((PetscObject)B)->type_name);CHKERRQ(ierr);
  ierr = MatSetUp(B);CHKERRQ(ierr);
  ierr = PetscNewLog(B,&chol);CHKERRQ(ierr);

  chol->Wrap    = MatWrapCholmod_seqsbaij;
  B->data       = chol;

  B->ops->getinfo                = MatGetInfo_External;
  B->ops->view                   = MatView_CHOLMOD;
  B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_CHOLMOD;
  B->ops->destroy                = MatDestroy_CHOLMOD;
  ierr                           = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_seqsbaij_cholmod);CHKERRQ(ierr);
  B->factortype                  = MAT_FACTOR_CHOLESKY;
  B->assembled                   = PETSC_TRUE; /* required by -ksp_view */
  B->preallocated                = PETSC_TRUE;

  ierr = CholmodStart(B);CHKERRQ(ierr);

  ierr = PetscFree(B->solvertype);CHKERRQ(ierr);
  ierr = PetscStrallocpy(MATSOLVERCHOLMOD,&B->solvertype);CHKERRQ(ierr);

  *F   = B;
  PetscFunctionReturn(0);
}

/*MC
     KSPFBCGSR - Implements a mathematically equivalent variant of FBiCGSTab.

   Options Database Keys:
.   see KSPSolve()

   Level: beginner

   Notes:
    Only allow right preconditioning

.seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPBICG, KSPFBCGSL, KSPSetPCSide()
M*/
PETSC_EXTERN PetscErrorCode KSPCreate_FBCGSR(KSP ksp)
{
  PetscErrorCode ierr;
  KSP_BCGS       *bcgs;

  PetscFunctionBegin;
  ierr = PetscNewLog(ksp,&bcgs);CHKERRQ(ierr);

  ksp->data                = bcgs;
  ksp->ops->setup          = KSPSetUp_FBCGSR;
  ksp->ops->solve          = KSPSolve_FBCGSR;
  ksp->ops->destroy        = KSPDestroy_BCGS;
  ksp->ops->reset          = KSPReset_BCGS;
  ksp->ops->buildsolution  = KSPBuildSolution_BCGS;
  ksp->ops->buildresidual  = KSPBuildResidualDefault;
  ksp->ops->setfromoptions = KSPSetFromOptions_BCGS;
  ksp->pc_side             = PC_RIGHT; /* set default PC side */

  ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);CHKERRQ(ierr);
  ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PETSC_EXTERN PetscErrorCode KSPCreate_CG(KSP ksp)
{
  PetscErrorCode ierr;
  KSP_CG         *cg;

  PetscFunctionBegin;
  ierr = PetscNewLog(ksp,&cg);CHKERRQ(ierr);
#if !defined(PETSC_USE_COMPLEX)
  cg->type = KSP_CG_SYMMETRIC;
#else
  cg->type = KSP_CG_HERMITIAN;
#endif
  ksp->data = (void*)cg;

  ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);CHKERRQ(ierr);
  ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);CHKERRQ(ierr);
  ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);CHKERRQ(ierr);

  /*
       Sets the functions that are associated with this data structure
       (in C++ this is the same as defining virtual functions)
  */
  ksp->ops->setup          = KSPSetUp_CG;
  ksp->ops->solve          = KSPSolve_CG;
  ksp->ops->destroy        = KSPDestroy_CG;
  ksp->ops->view           = KSPView_CG;
  ksp->ops->setfromoptions = KSPSetFromOptions_CG;
  ksp->ops->buildsolution  = KSPBuildSolutionDefault;
  ksp->ops->buildresidual  = KSPBuildResidualDefault;

  /*
      Attach the function KSPCGSetType_CG() to this object. The routine
      KSPCGSetType() checks for this attached function and calls it if it finds
      it. (Sort of like a dynamic member function that can be added at run time
  */
  ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPCGSetType_C",KSPCGSetType_CG);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPCGUseSingleReduction_C",KSPCGUseSingleReduction_CG);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PETSC_EXTERN PetscErrorCode DMCreate_ADDA(DM dm)
{
  PetscErrorCode ierr;
  DM_ADDA        *dd;

  PetscFunctionBegin;
  ierr     = PetscNewLog(dm,DM_ADDA,&dd);CHKERRQ(ierr);
  dm->data = (void*)dd;

  dm->ops->view                = DMView;
  dm->ops->createglobalvector  = DMCreateGlobalVector_ADDA;
  dm->ops->getcoloring         = DMCreateColoring_ADDA;
  dm->ops->creatematrix        = DMCreateMatrix_ADDA;
  dm->ops->createinterpolation = DMCreateInterpolation_ADDA;
  dm->ops->refine              = DMRefine_ADDA;
  dm->ops->coarsen             = DMCoarsen_ADDA;
  dm->ops->getinjection        = DMCreateInjection_ADDA;
  dm->ops->getaggregates       = DMCreateAggregates_ADDA;
  dm->ops->setup               = DMSetUp_ADDA;
  dm->ops->destroy             = DMDestroy_ADDA;
  PetscFunctionReturn(0);
}

EXTERN_C_END

EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "MatGetFactor_mpisbaij_pastix"
PetscErrorCode MatGetFactor_mpisbaij_pastix(Mat A,MatFactorType ftype,Mat *F)
{
  Mat            B;
  PetscErrorCode ierr;
  Mat_Pastix    *pastix;

  PetscFunctionBegin;
  if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with PaStiX LU, use AIJ matrix");

  /* Create the factorization matrix */
  ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr);
  ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr);
  ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr);
  ierr = MatSeqSBAIJSetPreallocation(B,1,0,PETSC_NULL);CHKERRQ(ierr);
  ierr = MatMPISBAIJSetPreallocation(B,1,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr);

  B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJPASTIX;
  B->ops->view                   = MatView_PaStiX;
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_pastix",MatFactorGetSolverPackage_pastix);CHKERRQ(ierr);
  B->factortype                  = MAT_FACTOR_CHOLESKY;

  ierr = PetscNewLog(B,Mat_Pastix,&pastix);CHKERRQ(ierr);
  pastix->CleanUpPastix             = PETSC_FALSE;
  pastix->isAIJ                     = PETSC_TRUE;
  pastix->scat_rhs                  = PETSC_NULL;
  pastix->scat_sol                  = PETSC_NULL;
  pastix->Destroy                   = B->ops->destroy;
  B->ops->destroy                   = MatDestroy_Pastix;
  B->spptr                          = (void*)pastix;

  *F = B;
  PetscFunctionReturn(0);
}

EXTERN_C_END

EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "MatGetFactor_seqaij_cholmod"
/* Almost a copy of MatGetFactor_seqsbaij_cholmod, yuck */
PetscErrorCode MatGetFactor_seqaij_cholmod(Mat A,MatFactorType ftype,Mat *F)
{
  Mat            B;
  Mat_CHOLMOD    *chol;
  PetscErrorCode ierr;
  PetscInt       m=A->rmap->n,n=A->cmap->n;

  PetscFunctionBegin;
  if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"CHOLMOD cannot do %s factorization with AIJ, only %s",
                                             MatFactorTypes[ftype],MatFactorTypes[MAT_FACTOR_CHOLESKY]);
  /* Create the factorization matrix F */
  ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr);
  ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr);
  ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr);
  ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr);
  ierr = PetscNewLog(B,Mat_CHOLMOD,&chol);CHKERRQ(ierr);
  chol->Wrap    = MatWrapCholmod_seqaij;
  chol->Destroy = MatDestroy_SeqAIJ;
  B->spptr      = chol;

  B->ops->view                   = MatView_CHOLMOD;
  B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_CHOLMOD;
  B->ops->destroy                = MatDestroy_CHOLMOD;
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_seqaij_cholmod",MatFactorGetSolverPackage_seqaij_cholmod);CHKERRQ(ierr);
  B->factortype   = MAT_FACTOR_CHOLESKY;
  B->assembled    = PETSC_TRUE;  /* required by -ksp_view */
  B->preallocated = PETSC_TRUE;

  ierr = CholmodStart(B);CHKERRQ(ierr);
  *F = B;
  PetscFunctionReturn(0);
}

/* ---------------------------------------------------------- */
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "TaoCreate_NM"
PetscErrorCode TaoCreate_NM(Tao tao)
{
  TAO_NelderMead *nm;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = PetscNewLog(tao,&nm);CHKERRQ(ierr);
  tao->data = (void*)nm;

  tao->ops->setup = TaoSetUp_NM;
  tao->ops->solve = TaoSolve_NM;
  tao->ops->view = TaoView_NM;
  tao->ops->setfromoptions = TaoSetFromOptions_NM;
  tao->ops->destroy = TaoDestroy_NM;

  tao->max_it = 2000;
  tao->max_funcs = 4000;
#if defined(PETSC_USE_REAL_SINGLE)
  tao->fatol = 1e-5;
  tao->frtol = 1e-5;
#else
  tao->fatol = 1e-8;
  tao->frtol = 1e-8;
#endif

  nm->simplex = 0;
  nm->lamda = 1;

  nm->mu_ic = -0.5;
  nm->mu_oc = 0.5;
  nm->mu_r = 1.0;
  nm->mu_e = 2.0;

  PetscFunctionReturn(0);
}

PETSC_EXTERN PetscErrorCode KSPCreate_Chebyshev(KSP ksp)
{
  PetscErrorCode ierr;
  KSP_Chebyshev  *chebyshevP;

  PetscFunctionBegin;
  ierr = PetscNewLog(ksp,&chebyshevP);CHKERRQ(ierr);

  ksp->data = (void*)chebyshevP;
  ierr      = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);CHKERRQ(ierr);
  ierr      = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);CHKERRQ(ierr);

  chebyshevP->emin = 0.;
  chebyshevP->emax = 0.;

  chebyshevP->tform[0] = 0.0;
  chebyshevP->tform[1] = 0.1;
  chebyshevP->tform[2] = 0;
  chebyshevP->tform[3] = 1.1;
  chebyshevP->eststeps = 10;
  chebyshevP->userandom = PETSC_FALSE;

  ksp->ops->setup          = KSPSetUp_Chebyshev;
  ksp->ops->solve          = KSPSolve_Chebyshev;
  ksp->ops->destroy        = KSPDestroy_Chebyshev;
  ksp->ops->buildsolution  = KSPBuildSolutionDefault;
  ksp->ops->buildresidual  = KSPBuildResidualDefault;
  ksp->ops->setfromoptions = KSPSetFromOptions_Chebyshev;
  ksp->ops->view           = KSPView_Chebyshev;
  ksp->ops->reset          = KSPReset_Chebyshev;

  ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetEigenvalues_C",KSPChebyshevSetEigenvalues_Chebyshev);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSet_C",KSPChebyshevEstEigSet_Chebyshev);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSetRandom_C",KSPChebyshevEstEigSetRandom_Chebyshev);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSetUseRandom_C",KSPChebyshevEstEigSetUseRandom_Chebyshev);CHKERRQ(ierr);
  ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigGetKSP_C",KSPChebyshevEstEigGetKSP_Chebyshev);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}