int TaoGPCGCreateLineSearch(TAO_SOLVER tao)
{
  int info;
  TAO_GPCGLINESEARCH *neP;

  TaoFunctionBegin;

  info = TaoNew(TAO_GPCGLINESEARCH,&neP);CHKERRQ(info);
  info = PetscLogObjectMemory(tao,sizeof(TAO_GPCGLINESEARCH)); CHKERRQ(info);
  neP->ftol		  = 0.05;
  neP->rtol		  = 0.0;
  neP->gtol		  = 0.0;
  neP->stepmin		  = 1.0e-20;
  neP->stepmax		  = 1.0e+20;
  neP->nfev		  = 0; 
  neP->bracket		  = 0; 
  neP->infoc              = 1;
  neP->maxfev		  = 30;
  neP->setupcalled        = 0;

  info = TaoSetLineSearch(tao,0,
			  TaoGPCGSetOptionsLineSearch,
			  TaoGPCGApplyLineSearch,
			  TaoGPCGViewLineSearch,
			  TaoGPCGDestroyLineSearch,
			  (void *) neP);CHKERRQ(info);

  TaoFunctionReturn(0);
}

PetscErrorCode    KSPSetUp_PIPEFCG(KSP ksp)
{
  PetscErrorCode ierr;
  KSP_PIPEFCG    *pipefcg;
  const PetscInt nworkstd = 5;

  PetscFunctionBegin;
  pipefcg = (KSP_PIPEFCG*)ksp->data;

  /* Allocate "standard" work vectors (not including the basis and transformed basis vectors) */
  ierr = KSPSetWorkVecs(ksp,nworkstd);CHKERRQ(ierr);

  /* Allocated space for pointers to additional work vectors
   note that mmax is the number of previous directions, so we add 1 for the current direction,
   and an extra 1 for the prealloc (which might be empty) */
  ierr = PetscMalloc4(pipefcg->mmax+1,&(pipefcg->Pvecs),pipefcg->mmax+1,&(pipefcg->pPvecs),pipefcg->mmax+1,&(pipefcg->Svecs),pipefcg->mmax+1,&(pipefcg->pSvecs));CHKERRQ(ierr);
  ierr = PetscMalloc4(pipefcg->mmax+1,&(pipefcg->Qvecs),pipefcg->mmax+1,&(pipefcg->pQvecs),pipefcg->mmax+1,&(pipefcg->ZETAvecs),pipefcg->mmax+1,&(pipefcg->pZETAvecs));CHKERRQ(ierr);
  ierr = PetscMalloc4(pipefcg->mmax+1,&(pipefcg->Pold),pipefcg->mmax+1,&(pipefcg->Sold),pipefcg->mmax+1,&(pipefcg->Qold),pipefcg->mmax+1,&(pipefcg->ZETAold));CHKERRQ(ierr);
  ierr = PetscMalloc1(pipefcg->mmax+1,&(pipefcg->chunksizes));CHKERRQ(ierr);
  ierr = PetscMalloc3(pipefcg->mmax+2,&(pipefcg->dots),pipefcg->mmax+1,&(pipefcg->etas),pipefcg->mmax+2,&(pipefcg->redux));CHKERRQ(ierr);

  /* If the requested number of preallocated vectors is greater than mmax reduce nprealloc */
  if(pipefcg->nprealloc > pipefcg->mmax+1){
    ierr = PetscInfo2(NULL,"Requested nprealloc=%d is greater than m_max+1=%d. Resetting nprealloc = m_max+1.\n",pipefcg->nprealloc, pipefcg->mmax+1);CHKERRQ(ierr);
  }

  /* Preallocate additional work vectors */
  ierr = KSPAllocateVectors_PIPEFCG(ksp,pipefcg->nprealloc,pipefcg->nprealloc);CHKERRQ(ierr);

  ierr = PetscLogObjectMemory((PetscObject)ksp,(pipefcg->mmax+1)*4*sizeof(Vec*)+(pipefcg->mmax+1)*4*sizeof(Vec**)+(pipefcg->mmax+1)*4*sizeof(Vec*)+
    (pipefcg->mmax+1)*sizeof(PetscInt)+(pipefcg->mmax+2)*sizeof(Vec*)+(pipefcg->mmax+2)*sizeof(PetscScalar)+(pipefcg->mmax+1)*sizeof(PetscReal));CHKERRQ(ierr);
  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);
}

PetscErrorCode BVAllocateWork_Private(BV bv,PetscInt s)
{
  PetscErrorCode ierr;

  PetscFunctionBegin;
  if (s>bv->lwork) {
    ierr = PetscFree(bv->work);CHKERRQ(ierr);
    ierr = PetscMalloc1(s,&bv->work);CHKERRQ(ierr);
    ierr = PetscLogObjectMemory((PetscObject)bv,(s-bv->lwork)*sizeof(PetscScalar));CHKERRQ(ierr);
    bv->lwork = s;
  }
  PetscFunctionReturn(0);
}

PetscErrorCode KSPSetUp_LCD(KSP ksp)
{
  KSP_LCD        *lcd = (KSP_LCD*)ksp->data;
  PetscErrorCode ierr;
  PetscInt       restart = lcd->restart;

  PetscFunctionBegin;
  /* get work vectors needed by LCD */
  ierr = KSPSetWorkVecs(ksp,2);CHKERRQ(ierr);

  ierr = VecDuplicateVecs(ksp->work[0],restart+1,&lcd->P);CHKERRQ(ierr);
  ierr = VecDuplicateVecs(ksp->work[0], restart + 1, &lcd->Q);CHKERRQ(ierr);
  ierr = PetscLogObjectMemory((PetscObject)ksp,2*(restart+2)*sizeof(Vec));CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*@
   SVDAllocateSolution - Allocate memory storage for common variables such
   as the singular values and the basis vectors.

   Collective on SVD

   Input Parameters:
+  svd   - eigensolver context
-  extra - number of additional positions, used for methods that require a
           working basis slightly larger than ncv

   Developers Notes:
   This is PETSC_EXTERN because it may be required by user plugin SVD
   implementations.

   This is called at setup after setting the value of ncv and the flag leftbasis.

   Level: developer
@*/
PetscErrorCode SVDAllocateSolution(SVD svd,PetscInt extra)
{
  PetscErrorCode ierr;
  PetscInt       oldsize,requested;
  Vec            tr,tl;

  PetscFunctionBegin;
  requested = svd->ncv + extra;

  /* oldsize is zero if this is the first time setup is called */
  ierr = BVGetSizes(svd->V,NULL,NULL,&oldsize);CHKERRQ(ierr);

  /* allocate sigma */
  if (requested != oldsize) {
    if (oldsize) {
      ierr = PetscFree3(svd->sigma,svd->perm,svd->errest);CHKERRQ(ierr);
    }
    ierr = PetscMalloc3(requested,&svd->sigma,requested,&svd->perm,requested,&svd->errest);CHKERRQ(ierr);
    ierr = PetscLogObjectMemory((PetscObject)svd,PetscMax(0,requested-oldsize)*(2*sizeof(PetscReal)+sizeof(PetscInt)));CHKERRQ(ierr);
  }
  /* allocate V */
  if (!svd->V) { ierr = SVDGetBV(svd,&svd->V,NULL);CHKERRQ(ierr); }
  if (!oldsize) {
    if (!((PetscObject)(svd->V))->type_name) {
      ierr = BVSetType(svd->V,BVSVEC);CHKERRQ(ierr);
    }
    ierr = SVDMatGetVecs(svd,&tr,NULL);CHKERRQ(ierr);
    ierr = BVSetSizesFromVec(svd->V,tr,requested);CHKERRQ(ierr);
    ierr = VecDestroy(&tr);CHKERRQ(ierr);
  } else {
    ierr = BVResize(svd->V,requested,PETSC_FALSE);CHKERRQ(ierr);
  }
  /* allocate U */
  if (svd->leftbasis) {
    if (!svd->U) { ierr = SVDGetBV(svd,NULL,&svd->U);CHKERRQ(ierr); }
    if (!oldsize) {
      if (!((PetscObject)(svd->U))->type_name) {
        ierr = BVSetType(svd->U,BVSVEC);CHKERRQ(ierr);
      }
      ierr = SVDMatGetVecs(svd,NULL,&tl);CHKERRQ(ierr);
      ierr = BVSetSizesFromVec(svd->U,tl,requested);CHKERRQ(ierr);
      ierr = VecDestroy(&tl);CHKERRQ(ierr);
    } else {
      ierr = BVResize(svd->U,requested,PETSC_FALSE);CHKERRQ(ierr);
    }
  }
  PetscFunctionReturn(0);
}

/*@C
   PetscDrawHGCreate - Creates a histogram data structure.

   Collective over PetscDraw

   Input Parameters:
+  draw  - The window where the graph will be made
-  bins - The number of bins to use

   Output Parameters:
.  hist - The histogram context

   Level: intermediate

   Concepts: histogram^creating

.seealso: PetscDrawHGDestroy()

@*/
PetscErrorCode  PetscDrawHGCreate(PetscDraw draw, int bins, PetscDrawHG *hist)
{
  PetscBool      isnull;
  PetscDrawHG    h;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(draw, PETSC_DRAW_CLASSID,1);
  PetscValidLogicalCollectiveInt(draw,bins,2);
  PetscValidPointer(hist,3);

  ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr);
  if (isnull) {*hist = NULL; PetscFunctionReturn(0);}

  ierr = PetscHeaderCreate(h, PETSC_DRAWHG_CLASSID,  "PetscDrawHG", "Histogram", "Draw", PetscObjectComm((PetscObject)draw), PetscDrawHGDestroy, NULL);CHKERRQ(ierr);
  ierr = PetscLogObjectParent((PetscObject)draw,(PetscObject)h);CHKERRQ(ierr);

  ierr = PetscObjectReference((PetscObject)draw);CHKERRQ(ierr);
  h->win = draw;

  h->view        = NULL;
  h->destroy     = NULL;
  h->color       = PETSC_DRAW_GREEN;
  h->xmin        = PETSC_MAX_REAL;
  h->xmax        = PETSC_MIN_REAL;
  h->ymin        = 0.;
  h->ymax        = 1.;
  h->numBins     = bins;
  h->maxBins     = bins;

  ierr = PetscMalloc1(h->maxBins,&h->bins);CHKERRQ(ierr);

  h->numValues   = 0;
  h->maxValues   = CHUNKSIZE;
  h->calcStats   = PETSC_FALSE;
  h->integerBins = PETSC_FALSE;

  ierr = PetscMalloc1(h->maxValues,&h->values);CHKERRQ(ierr);
  ierr = PetscLogObjectMemory((PetscObject)h,(h->maxBins + h->maxValues)*sizeof(PetscReal));CHKERRQ(ierr);

  ierr = PetscDrawAxisCreate(draw,&h->axis);CHKERRQ(ierr);
  ierr = PetscLogObjectParent((PetscObject)h,(PetscObject)h->axis);CHKERRQ(ierr);

  *hist = h;
  PetscFunctionReturn(0);
}

EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "KSPCreate_BSSCR"
PetscErrorCode PETSCKSP_DLLEXPORT KSPCreate_BSSCR(KSP ksp)
{
    KSP_BSSCR  *bsscr;
    PetscErrorCode ierr;
    PetscFunctionBegin;
    ierr = Stg_PetscNew(KSP_BSSCR,&bsscr);CHKERRQ(ierr);
    ierr = PetscLogObjectMemory((PetscObject)ksp,sizeof(KSP_BSSCR));CHKERRQ(ierr);
    //ierr = PetscNewLog(ksp,KSP_BSSCR,&bsscr);CHKERRQ(ierr);
    ksp->data                              = (void*)bsscr;

    #if ( (PETSC_VERSION_MAJOR >= 3) && (PETSC_VERSION_MINOR >= 2 ) )
    ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,0);CHKERRQ(ierr);
    ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,1);CHKERRQ(ierr);
    ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,0);CHKERRQ(ierr);
    ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);CHKERRQ(ierr);
    #endif
    /*
       Sets the functions that are associated with this data structure
       (in C++ this is the same as defining virtual functions)
    */
    ksp->ops->setup                = KSPSetUp_BSSCR;
    ksp->ops->solve                = KSPSolve_BSSCR;
    ksp->ops->destroy              = KSPDestroy_BSSCR;
    ksp->ops->view                 = KSPView_BSSCR;
    ksp->ops->setfromoptions       = KSPSetFromOptions_BSSCR;
    ksp->ops->buildsolution        = KSPDefaultBuildSolution;
    ksp->ops->buildresidual        = KSPDefaultBuildResidual;

//    bsscr->k2type=K2_GMG;
    bsscr->k2type      = 0;
    bsscr->do_scaling  = PETSC_FALSE;
    bsscr->scaled      = PETSC_FALSE;
    bsscr->K2built     = PETSC_FALSE;
    bsscr->check_cb_pressureNS     = PETSC_FALSE;/* checker board nullspace */
    bsscr->check_const_pressureNS  = PETSC_FALSE;/* constant nullspace */
    bsscr->check_pressureNS        = PETSC_FALSE;/* is true if either of above two are true */
    bsscr->t           = NULL;/* null space vectors for pressure */
    bsscr->v           = NULL;
    bsscr->nstol       = 1e-7;/* null space detection tolerance */
    bsscr->uStar       = NULL;
    bsscr->been_here   = 0;
    PetscFunctionReturn(0);
}

/*@
   PetscDrawHGSetNumberBins - Change the number of bins that are to be drawn.

   Not Collective (ignored except on processor 0 of PetscDrawHG)

   Input Parameter:
+  hist - The histogram context.
-  bins  - The number of bins.

   Level: intermediate

   Concepts: histogram^setting number of bins

@*/
PetscErrorCode  PetscDrawHGSetNumberBins(PetscDrawHG hist, int bins)
{
  PetscErrorCode ierr;

  PetscFunctionBegin;
  if (!hist) PetscFunctionReturn(0);
  PetscValidHeaderSpecific(hist,PETSC_DRAWHG_CLASSID,1);
  PetscValidLogicalCollectiveInt(hist,bins,2);

  if (hist->maxBins < bins) {
    ierr = PetscFree(hist->bins);CHKERRQ(ierr);
    ierr = PetscMalloc1(bins, &hist->bins);CHKERRQ(ierr);
    ierr = PetscLogObjectMemory((PetscObject)hist, (bins - hist->maxBins) * sizeof(PetscReal));CHKERRQ(ierr);
    hist->maxBins = bins;
  }
  hist->numBins = bins;
  PetscFunctionReturn(0);
}

PETSC_EXTERN PetscErrorCode CharacteristicCreate_DA(Characteristic c)
{
  Characteristic_DA *da;
  PetscErrorCode    ierr;

  PetscFunctionBegin;
  ierr    = PetscNew(&da);CHKERRQ(ierr);
  ierr    = PetscMemzero(da, sizeof(Characteristic_DA));CHKERRQ(ierr);
  ierr    = PetscLogObjectMemory((PetscObject)c, sizeof(Characteristic_DA));CHKERRQ(ierr);
  c->data = (void*) da;

  c->ops->setup   = CharacteristicSetUp_DA;
  c->ops->destroy = CharacteristicDestroy_DA;
  c->ops->view    = CharacteristicView_DA;

  da->dummy = 0;
  PetscFunctionReturn(0);
}

PetscErrorCode BSSCR_PCCreate_GtKG( PC pc )
{
	PC_GtKG         pc_data;
	PetscErrorCode  ierr;
	
	/* create memory for ctx */
	ierr = Stg_PetscNew( _PC_GtKG,&pc_data);CHKERRQ(ierr);
	
	/* init ctx */
	pc_data->K   = PETSC_NULL;
	pc_data->G   = PETSC_NULL;
	pc_data->M   = PETSC_NULL;
	pc_data->GtG = PETSC_NULL;
	pc_data->form_GtG          = PETSC_TRUE;
	pc_data->ksp               = PETSC_NULL;
	pc_data->monitor_activated = PETSC_FALSE;
	pc_data->monitor_rhs_consistency = PETSC_FALSE;
	
	pc_data->s = PETSC_NULL;
	pc_data->t = PETSC_NULL;
	pc_data->X = PETSC_NULL;
	pc_data->inv_diag_M = PETSC_NULL;
	
	/* create internals */
	KSPCreate( ((PetscObject)pc)->comm, &pc_data->ksp );
	
	
	/* set ctx onto pc */
	pc->data  = (void*)pc_data;
	
	ierr = PetscLogObjectMemory(pc,sizeof(_PC_GtKG));CHKERRQ(ierr);
	
	/* define operations */
	pc->ops->setup   = BSSCR_PCSetUp_GtKG;
	pc->ops->view    = BSSCR_PCView_GtKG;
	pc->ops->destroy = BSSCR_PCDestroy_GtKG;
	pc->ops->setfromoptions = BSSCR_PCSetFromOptions_GtKG;
	
	pc->ops->apply          = BSSCR_PCApply_GtKG;
	pc->ops->applytranspose = BSSCR_PCApplyTranspose_GtKG;
	
	
	PetscFunctionReturn(0);
}

PetscErrorCode BSSCR_PCCreate_ScGtKG( PC pc )
{
	PC_SC_GtKG      pc_data;
	PetscErrorCode  ierr;
	
	/* create memory for ctx */
	ierr = Stg_PetscNew( _PC_SC_GtKG,&pc_data);CHKERRQ(ierr);
	
	/* init ctx */
	pc_data->F   = PETSC_NULL;
	pc_data->Bt  = PETSC_NULL;
	pc_data->B   = PETSC_NULL;
	pc_data->BBt_has_cnst_nullspace = PETSC_FALSE;
	pc_data->ksp_BBt   = PETSC_NULL;
	pc_data->monitor_activated = PETSC_FALSE;
	
	pc_data->X1 = PETSC_NULL;
	pc_data->X2 = PETSC_NULL;
	pc_data->Y1 = PETSC_NULL;
	pc_data->Y2 = PETSC_NULL;
	
	pc_data->s = PETSC_NULL;
	pc_data->t = PETSC_NULL;
	pc_data->X = PETSC_NULL;
	
	
	
	/* set ctx onto pc */
	pc->data  = (void*)pc_data;
	
	ierr = PetscLogObjectMemory(pc,sizeof(_PC_SC_GtKG));CHKERRQ(ierr);
	
	/* define operations */
	pc->ops->setup   = BSSCR_PCSetUp_ScGtKG;
	pc->ops->view    = BSSCR_PCView_ScGtKG;
	pc->ops->destroy = BSSCR_PCDestroy_ScGtKG;
	pc->ops->setfromoptions = BSSCR_PCSetFromOptions_ScGtKG;
	
	pc->ops->apply          = BSSCR_PCApply_ScGtKG;
	pc->ops->applytranspose = BSSCR_PCApplyTranspose_ScGtKG;
	
	
	PetscFunctionReturn(0);
}

PetscErrorCode  KSPFischerGuessCreate_Method2(KSP ksp,int  maxl,KSPFischerGuess_Method2 **ITG)
{
  KSPFischerGuess_Method2 *itg;
  PetscErrorCode          ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
  ierr = PetscMalloc(sizeof(KSPFischerGuess_Method2),&itg);CHKERRQ(ierr);
  ierr = PetscMalloc(maxl * sizeof(PetscScalar),&itg->alpha);CHKERRQ(ierr);
  ierr = PetscLogObjectMemory(ksp,sizeof(KSPFischerGuess_Method2) + maxl*sizeof(PetscScalar));CHKERRQ(ierr);
  ierr = KSPGetVecs(ksp,maxl,&itg->xtilde,0,PETSC_NULL);CHKERRQ(ierr);
  ierr = PetscLogObjectParents(ksp,maxl,itg->xtilde);CHKERRQ(ierr);
  ierr = VecDuplicate(itg->xtilde[0],&itg->Ax);CHKERRQ(ierr);
  ierr = PetscLogObjectParent(ksp,itg->Ax);CHKERRQ(ierr);
  ierr = VecDuplicate(itg->xtilde[0],&itg->guess);CHKERRQ(ierr);
  ierr = PetscLogObjectParent(ksp,itg->guess);CHKERRQ(ierr);
  *ITG = itg;
  PetscFunctionReturn(0);
}

/*
    VecCreate_MPI_Private - Basic create routine called by VecCreate_MPI() (i.e. VecCreateMPI()),
    VecCreateMPIWithArray(), VecCreate_Shared() (i.e. VecCreateShared()), VecCreateGhost(),
    VecDuplicate_MPI(), VecCreateGhostWithArray(), VecDuplicate_MPI(), and VecDuplicate_Shared()

    If alloc is true and array is NULL then this routine allocates the space, otherwise
    no space is allocated.
*/
PetscErrorCode VecCreate_MPI_Private(Vec v,PetscBool alloc,PetscInt nghost,const PetscScalar array[])
{
  Vec_MPI        *s;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr           = PetscNewLog(v,&s);CHKERRQ(ierr);
  v->data        = (void*)s;
  ierr           = PetscMemcpy(v->ops,&DvOps,sizeof(DvOps));CHKERRQ(ierr);
  s->nghost      = nghost;
  v->petscnative = PETSC_TRUE;

  ierr = PetscLayoutSetUp(v->map);CHKERRQ(ierr);

  s->array           = (PetscScalar*)array;
  s->array_allocated = 0;
  if (alloc && !array) {
    PetscInt n = v->map->n+nghost;
    ierr               = PetscMalloc1(n,&s->array);CHKERRQ(ierr);
    ierr               = PetscLogObjectMemory((PetscObject)v,n*sizeof(PetscScalar));CHKERRQ(ierr);
    ierr               = PetscMemzero(s->array,n*sizeof(PetscScalar));CHKERRQ(ierr);
    s->array_allocated = s->array;
  }

  /* By default parallel vectors do not have local representation */
  s->localrep    = 0;
  s->localupdate = 0;

  v->stash.insertmode = NOT_SET_VALUES;
  v->bstash.insertmode = NOT_SET_VALUES;
  /* create the stashes. The block-size for bstash is set later when
     VecSetValuesBlocked is called.
  */
  ierr = VecStashCreate_Private(PetscObjectComm((PetscObject)v),1,&v->stash);CHKERRQ(ierr);
  ierr = VecStashCreate_Private(PetscObjectComm((PetscObject)v),PetscAbs(v->map->bs),&v->bstash);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
  ierr = PetscObjectChangeTypeName((PetscObject)v,VECMPI);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*@
   NEPAllocateSolution - Allocate memory storage for common variables such
   as eigenvalues and eigenvectors.

   Collective on NEP

   Input Parameters:
+  nep   - eigensolver context
-  extra - number of additional positions, used for methods that require a
           working basis slightly larger than ncv

   Developers Note:
   This is PETSC_EXTERN because it may be required by user plugin NEP
   implementations.

   Level: developer
@*/
PetscErrorCode NEPAllocateSolution(NEP nep,PetscInt extra)
{
  PetscErrorCode ierr;
  PetscInt       oldsize,newc,requested;
  PetscLogDouble cnt;
  Mat            T;
  Vec            t;

  PetscFunctionBegin;
  requested = nep->ncv + extra;

  /* oldsize is zero if this is the first time setup is called */
  ierr = BVGetSizes(nep->V,NULL,NULL,&oldsize);CHKERRQ(ierr);
  newc = PetscMax(0,requested-oldsize);

  /* allocate space for eigenvalues and friends */
  if (requested != oldsize) {
    if (oldsize) {
      ierr = PetscFree4(nep->eigr,nep->eigi,nep->errest,nep->perm);CHKERRQ(ierr);
    }
    ierr = PetscMalloc4(requested,&nep->eigr,requested,&nep->eigi,requested,&nep->errest,requested,&nep->perm);CHKERRQ(ierr);
    cnt = newc*sizeof(PetscScalar) + newc*sizeof(PetscReal) + newc*sizeof(PetscInt);
    ierr = PetscLogObjectMemory((PetscObject)nep,cnt);CHKERRQ(ierr);
  }

  /* allocate V */
  if (!nep->V) { ierr = NEPGetBV(nep,&nep->V);CHKERRQ(ierr); }
  if (!oldsize) {
    if (!((PetscObject)(nep->V))->type_name) {
      ierr = BVSetType(nep->V,BVSVEC);CHKERRQ(ierr);
    }
    if (nep->split) T = nep->A[0];
    else {
      ierr = NEPGetFunction(nep,&T,NULL,NULL,NULL);CHKERRQ(ierr);
    }
    ierr = MatGetVecs(T,&t,NULL);CHKERRQ(ierr);
    ierr = BVSetSizesFromVec(nep->V,t,requested);CHKERRQ(ierr);
    ierr = VecDestroy(&t);CHKERRQ(ierr);
  } else {
    ierr = BVResize(nep->V,requested,PETSC_FALSE);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

PetscErrorCode  KSPFischerGuessCreate_Method1(KSP ksp,int maxl,KSPFischerGuess_Method1 **ITG)
{
  KSPFischerGuess_Method1 *itg;
  PetscErrorCode          ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
  ierr = PetscNew(&itg);CHKERRQ(ierr);
  ierr = PetscMalloc1(maxl,&itg->alpha);CHKERRQ(ierr);
  ierr = PetscLogObjectMemory((PetscObject)ksp,sizeof(KSPFischerGuess_Method1) + maxl*sizeof(PetscScalar));CHKERRQ(ierr);
  ierr = KSPCreateVecs(ksp,maxl,&itg->xtilde,0,NULL);CHKERRQ(ierr);
  ierr = PetscLogObjectParents(ksp,maxl,itg->xtilde);CHKERRQ(ierr);
  ierr = KSPCreateVecs(ksp,maxl,&itg->btilde,0,NULL);CHKERRQ(ierr);
  ierr = PetscLogObjectParents(ksp,maxl,itg->btilde);CHKERRQ(ierr);
  ierr = VecDuplicate(itg->xtilde[0],&itg->guess);CHKERRQ(ierr);
  ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)itg->guess);CHKERRQ(ierr);
  *ITG = itg;
  PetscFunctionReturn(0);
}

/*------------------------------------------------------------*/
EXTERN_C_BEGIN
#undef __FUNCT__  
#define __FUNCT__ "TaoCreate_GPCG"
int TaoCreate_GPCG(TAO_SOLVER tao)
{
  TAO_GPCG *gpcg;
  int      info;

  TaoFunctionBegin;

  info = TaoNew(TAO_GPCG,&gpcg); CHKERRQ(info);
  info = PetscLogObjectMemory(tao,sizeof(TAO_GPCG)); CHKERRQ(info);

  info=TaoSetTaoSolveRoutine(tao,TaoSolve_GPCG,(void*)gpcg); CHKERRQ(info);
  info=TaoSetTaoSetUpDownRoutines(tao,TaoSetUp_GPCG,TaoSetDown_GPCG); CHKERRQ(info);
  info=TaoSetTaoOptionsRoutine(tao,TaoSetFromOptions_GPCG); CHKERRQ(info);
  info=TaoSetTaoViewRoutine(tao,TaoView_GPCG); CHKERRQ(info);
  info=TaoSetTaoDualVariablesRoutine(tao,TaoGetDualVariables_GPCG); CHKERRQ(info);

  info = TaoSetMaximumIterates(tao,500); CHKERRQ(info);
  info = TaoSetMaximumFunctionEvaluations(tao,100000); CHKERRQ(info);
  info = TaoSetTolerances(tao,1e-12,1e-12,0,0); CHKERRQ(info);

  /* Initialize pointers and variables */
  gpcg->n=0;
  gpcg->maxgpits = 8;
  gpcg->pg_ftol = 0.1;

  gpcg->gp_iterates=0; /* Cumulative number */
  gpcg->total_gp_its = 0;
 
  /* Initialize pointers and variables */
  gpcg->n_bind=0;
  gpcg->n_free = 0;
  gpcg->n_upper=0;
  gpcg->n_lower=0;

  //  info = TaoCreateProjectedLineSearch(tao); CHKERRQ(info);
  info = TaoGPCGCreateLineSearch(tao); CHKERRQ(info);

  TaoFunctionReturn(0);
}

/*@
   STSetUp - Prepares for the use of a spectral transformation.

   Collective on ST

   Input Parameter:
.  st - the spectral transformation context

   Level: advanced

.seealso: STCreate(), STApply(), STDestroy()
@*/
PetscErrorCode STSetUp(ST st)
{
  PetscInt       i,n,k;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(st,ST_CLASSID,1);
  STCheckMatrices(st,1);
  if (st->setupcalled) PetscFunctionReturn(0);
  ierr = PetscInfo(st,"Setting up new ST\n");CHKERRQ(ierr);
  ierr = PetscLogEventBegin(ST_SetUp,st,0,0,0);CHKERRQ(ierr);
  if (!((PetscObject)st)->type_name) {
    ierr = STSetType(st,STSHIFT);CHKERRQ(ierr);
  }
  if (!st->T) {
    ierr = PetscMalloc(PetscMax(2,st->nmat)*sizeof(Mat),&st->T);CHKERRQ(ierr);
    ierr = PetscLogObjectMemory((PetscObject)st,PetscMax(2,st->nmat)*sizeof(Mat));CHKERRQ(ierr);
    for (i=0;i<PetscMax(2,st->nmat);i++) st->T[i] = NULL;
  } else {
    for (i=0;i<PetscMax(2,st->nmat);i++) {
      ierr = MatDestroy(&st->T[i]);CHKERRQ(ierr);
    }
  }
  ierr = MatDestroy(&st->P);CHKERRQ(ierr);
  if (!st->w) {
    ierr = MatGetVecs(st->A[0],&st->w,NULL);CHKERRQ(ierr);
    ierr = PetscLogObjectParent((PetscObject)st,(PetscObject)st->w);CHKERRQ(ierr);
  }
  if (st->D) {
    ierr = MatGetLocalSize(st->A[0],NULL,&n);CHKERRQ(ierr);
    ierr = VecGetLocalSize(st->D,&k);CHKERRQ(ierr);
    if (n != k) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Balance matrix has wrong dimension %D (should be %D)",k,n);
    if (!st->wb) {
      ierr = VecDuplicate(st->D,&st->wb);CHKERRQ(ierr);
      ierr = PetscLogObjectParent((PetscObject)st,(PetscObject)st->wb);CHKERRQ(ierr);
    }
  }
  if (st->ops->setup) { ierr = (*st->ops->setup)(st);CHKERRQ(ierr); }
  st->setupcalled = 1;
  ierr = PetscLogEventEnd(ST_SetUp,st,0,0,0);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/*@C
   PetscDrawLGAddPoints - Adds several points to each of the line graphs.
   The new points must have an X coordinate larger than the old points.

   Not Collective, but ignored by all processors except processor 0 in PetscDrawLG

   Input Parameters:
+  lg - the LineGraph data structure
.  xx,yy - points to two arrays of pointers that point to arrays 
           containing the new x and y points for each curve.
-  n - number of points being added

   Level: intermediate


   Concepts: line graph^adding points

.seealso: PetscDrawLGAddPoint()
@*/
PetscErrorCode  PetscDrawLGAddPoints(PetscDrawLG lg,int n,PetscReal **xx,PetscReal **yy)
{
  PetscErrorCode ierr;
  int            i,j,k;
  PetscReal      *x,*y;

  PetscFunctionBegin;
  if (lg && ((PetscObject)lg)->classid == PETSC_DRAW_CLASSID) PetscFunctionReturn(0);
  PetscValidHeaderSpecific(lg,PETSC_DRAWLG_CLASSID,1);
  if (lg->loc+n*lg->dim >= lg->len) { /* allocate more space */
    PetscReal *tmpx,*tmpy;
    int    chunk = CHUNCKSIZE;

    if (n > chunk) chunk = n;
    ierr = PetscMalloc2(lg->len+lg->dim*chunk,PetscReal,&tmpx,lg->len+lg->dim*chunk,PetscReal,&tmpy);CHKERRQ(ierr);
    ierr = PetscLogObjectMemory(lg,2*lg->dim*chunk*sizeof(PetscReal));CHKERRQ(ierr);
    ierr = PetscMemcpy(tmpx,lg->x,lg->len*sizeof(PetscReal));CHKERRQ(ierr);
    ierr = PetscMemcpy(tmpy,lg->y,lg->len*sizeof(PetscReal));CHKERRQ(ierr);
    ierr = PetscFree2(lg->x,lg->y);CHKERRQ(ierr);
    lg->x = tmpx; 
    lg->y = tmpy;
    lg->len += lg->dim*chunk;
  }
  for (j=0; j<lg->dim; j++) {
    x = xx[j]; y = yy[j];
    k = lg->loc + j;
    for (i=0; i<n; i++) {
      if (x[i] > lg->xmax) lg->xmax = x[i]; 
      if (x[i] < lg->xmin) lg->xmin = x[i];
      if (y[i] > lg->ymax) lg->ymax = y[i]; 
      if (y[i] < lg->ymin) lg->ymin = y[i];

      lg->x[k]   = x[i];
      lg->y[k] = y[i];
      k += lg->dim;
    }
  }
  lg->loc   += n*lg->dim;
  lg->nopts += n;
  PetscFunctionReturn(0);
}

EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "CharacteristicCreate_DA"
PetscErrorCode CharacteristicCreate_DA(Characteristic c)
{
  Characteristic_DA *da;
  PetscErrorCode     ierr;

  PetscFunctionBegin;
  ierr = PetscNew(Characteristic_DA, &da);CHKERRQ(ierr);
  ierr = PetscMemzero(da, sizeof(Characteristic_DA));CHKERRQ(ierr);
  ierr = PetscLogObjectMemory(c, sizeof(Characteristic_DA));CHKERRQ(ierr);
  c->data = (void *) da;

  c->ops->setup   = CharacteristicSetUp_DA;
  c->ops->destroy = CharacteristicDestroy_DA;
  c->ops->view    = CharacteristicView_DA;

  da->dummy = 0;
  PetscFunctionReturn(0);
}