static PetscErrorCode PCApplyRichardson_MG(PC pc,Vec b,Vec x,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool zeroguess,PetscInt *outits,PCRichardsonConvergedReason *reason)
{
  PC_MG          *mg        = (PC_MG*)pc->data;
  PC_MG_Levels   **mglevels = mg->levels;
  PetscErrorCode ierr;
  PetscInt       levels = mglevels[0]->levels,i;

  PetscFunctionBegin;
  /* When the DM is supplying the matrix then it will not exist until here */
  for (i=0; i<levels; i++) {
    if (!mglevels[i]->A) {
      ierr = KSPGetOperators(mglevels[i]->smoothu,&mglevels[i]->A,NULL);CHKERRQ(ierr);
      ierr = PetscObjectReference((PetscObject)mglevels[i]->A);CHKERRQ(ierr);
    }
  }
  mglevels[levels-1]->b = b;
  mglevels[levels-1]->x = x;

  mg->rtol   = rtol;
  mg->abstol = abstol;
  mg->dtol   = dtol;
  if (rtol) {
    /* compute initial residual norm for relative convergence test */
    PetscReal rnorm;
    if (zeroguess) {
      ierr = VecNorm(b,NORM_2,&rnorm);CHKERRQ(ierr);
    } else {
      ierr = (*mglevels[levels-1]->residual)(mglevels[levels-1]->A,b,x,w);CHKERRQ(ierr);
      ierr = VecNorm(w,NORM_2,&rnorm);CHKERRQ(ierr);
    }
    mg->ttol = PetscMax(rtol*rnorm,abstol);
  } else if (abstol) mg->ttol = abstol;
  else mg->ttol = 0.0;

  /* since smoother is applied to full system, not just residual we need to make sure that smoothers don't
     stop prematurely due to small residual */
  for (i=1; i<levels; i++) {
    ierr = KSPSetTolerances(mglevels[i]->smoothu,0,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);
    if (mglevels[i]->smoothu != mglevels[i]->smoothd) {
      ierr = KSPSetTolerances(mglevels[i]->smoothd,0,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);
    }
  }

  *reason = (PCRichardsonConvergedReason)0;
  for (i=0; i<its; i++) {
    ierr = PCMGMCycle_Private(pc,mglevels+levels-1,reason);CHKERRQ(ierr);
    if (*reason) break;
  }
  if (!*reason) *reason = PCRICHARDSON_CONVERGED_ITS;
  *outits = i;
  PetscFunctionReturn(0);
}

/*@C
   KSPMonitorDynamicTolerance - Recompute the inner tolerance in every
   outer iteration in an adaptive way.

   Collective on KSP

   Input Parameters:
+  ksp   - iterative context
.  n     - iteration number (not used)
.  fnorm - the current residual norm
.  dummy - some context as a C struct. fields:
             coef: a scaling coefficient. default 1.0. can be passed through
                   -sub_ksp_dynamic_tolerance_param
             bnrm: norm of the right-hand side. store it to avoid repeated calculation

   Notes:
   This may be useful for a flexibly preconditioner Krylov method to
   control the accuracy of the inner solves needed to gaurantee the
   convergence of the outer iterations.

   Level: advanced

.keywords: KSP, inner tolerance

.seealso: KSPMonitorDynamicToleranceDestroy()
@*/
PetscErrorCode KSPMonitorDynamicTolerance(KSP ksp,PetscInt its,PetscReal fnorm,void *dummy)
{
    PetscErrorCode ierr;
    PC             pc;
    PetscReal      outer_rtol, outer_abstol, outer_dtol, inner_rtol;
    PetscInt       outer_maxits,nksp,first,i;
    KSPDynTolCtx   *scale   = (KSPDynTolCtx*)dummy;
    KSP            kspinner = NULL, *subksp = NULL;

    PetscFunctionBegin;
    ierr = KSPGetPC(ksp, &pc);
    CHKERRQ(ierr);

    /* compute inner_rtol */
    if (scale->bnrm < 0.0) {
        Vec b;
        ierr = KSPGetRhs(ksp, &b);
        CHKERRQ(ierr);
        ierr = VecNorm(b, NORM_2, &(scale->bnrm));
        CHKERRQ(ierr);
    }
    ierr       = KSPGetTolerances(ksp, &outer_rtol, &outer_abstol, &outer_dtol, &outer_maxits);
    CHKERRQ(ierr);
    inner_rtol = PetscMin(scale->coef * scale->bnrm * outer_rtol / fnorm, 0.999);
    /*ierr = PetscPrintf(PETSC_COMM_WORLD, "        Inner rtol = %g\n", (double)inner_rtol);CHKERRQ(ierr);*/

    /* if pc is ksp */
    ierr = PCKSPGetKSP(pc, &kspinner);
    CHKERRQ(ierr);
    if (kspinner) {
        ierr = KSPSetTolerances(kspinner, inner_rtol, outer_abstol, outer_dtol, outer_maxits);
        CHKERRQ(ierr);
        PetscFunctionReturn(0);
    }

    /* if pc is bjacobi */
    ierr = PCBJacobiGetSubKSP(pc, &nksp, &first, &subksp);
    CHKERRQ(ierr);
    if (subksp) {
        for (i=0; i<nksp; i++) {
            ierr = KSPSetTolerances(subksp[i], inner_rtol, outer_abstol, outer_dtol, outer_maxits);
            CHKERRQ(ierr);
        }
        PetscFunctionReturn(0);
    }

    /* todo: dynamic tolerance may apply to other types of pc too */
    PetscFunctionReturn(0);
}

MGSolver_PETScData* MultigridSolver_CreateOuterSolver( MultigridSolver* self, Mat matrix ) {
	//MatrixSolver*	outerSolver;
	MGSolver_PETScData* 	outerSolver = malloc( sizeof(MGSolver_PETScData) );
	PC			pc;

	/*
	outerSolver = (MatrixSolver*)PETScMatrixSolver_New( "" );
	PETScMatrixSolver_SetKSPType( outerSolver, PETScMatrixSolver_KSPType_Richardson );
	PETScMatrixSolver_SetPCType( outerSolver, PETScMatrixSolver_PCType_SOR );
	PETScMatrixSolver_SetMatrix( outerSolver, matrix );
	PETScMatrixSolver_SetMaxIterations( outerSolver, 3 );
	PETScMatrixSolver_SetUseInitialSolution( outerSolver, True );
	PETScMatrixSolver_SetNormType( outerSolver, PETScMatrixSolver_NormType_Preconditioned );
	*/
	KSPCreate( MPI_COMM_WORLD, &outerSolver->ksp );
	KSPSetType( outerSolver->ksp, KSPRICHARDSON );
	KSPGetPC( outerSolver->ksp, &pc );
	PCSetType( pc, PCSOR );
	if( outerSolver->matrix != PETSC_NULL )
		Stg_MatDestroy(&outerSolver->matrix );
	outerSolver->matrix = matrix;
	Stg_KSPSetOperators( outerSolver->ksp, matrix, matrix, DIFFERENT_NONZERO_PATTERN );
	KSPSetTolerances( outerSolver->ksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, (PetscInt)3 );
	KSPSetInitialGuessNonzero( outerSolver->ksp, (PetscTruth)True );
	KSPSetNormType( outerSolver->ksp, MultigridSolver_NormType_Preconditioned );

	return outerSolver;
}

int TaoLinearSolverPetsc::SetTolerances(double rtol, double atol,
					double dtol, int maxits) {
  int info;
  PetscFunctionBegin;
  info = KSPSetTolerances(ksp, rtol, atol, dtol, maxits); CHKERRQ(info);
  PetscFunctionReturn(0);
}

/*@C
   TSMonitorSPEigCtxCreate - Creates a context for use with TS to monitor the eigenvalues of the linearized operator

   Collective on TS

   Input Parameters:
+  host - the X display to open, or null for the local machine
.  label - the title to put in the title bar
.  x, y - the screen coordinates of the upper left coordinate of the window
.  m, n - the screen width and height in pixels
-  howoften - if positive then determines the frequency of the plotting, if -1 then only at the final time

   Output Parameter:
.  ctx - the context

   Options Database Key:
.  -ts_monitor_sp_eig - plot egienvalues of linearized right hand side

   Notes:
   Use TSMonitorSPEigCtxDestroy() to destroy.

   Currently only works if the Jacobian is provided explicitly.

   Currently only works for ODEs u_t - F(t,u) = 0; that is with no mass matrix.

   Level: intermediate

.keywords: TS, monitor, line graph, residual, seealso

.seealso: TSMonitorSPEigTimeStep(), TSMonitorSet(), TSMonitorLGSolution(), TSMonitorLGError()

@*/
PetscErrorCode  TSMonitorSPEigCtxCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscInt howoften,TSMonitorSPEigCtx *ctx)
{
  PetscDraw      win;
  PetscErrorCode ierr;
  PC             pc;

  PetscFunctionBegin;
  ierr = PetscNew(ctx);CHKERRQ(ierr);
  ierr = PetscRandomCreate(comm,&(*ctx)->rand);CHKERRQ(ierr);
  ierr = PetscRandomSetFromOptions((*ctx)->rand);CHKERRQ(ierr);
  ierr = PetscDrawCreate(comm,host,label,x,y,m,n,&win);CHKERRQ(ierr);
  ierr = PetscDrawSetFromOptions(win);CHKERRQ(ierr);
  ierr = PetscDrawSPCreate(win,1,&(*ctx)->drawsp);CHKERRQ(ierr);
  ierr = KSPCreate(comm,&(*ctx)->ksp);CHKERRQ(ierr);
  ierr = KSPSetOptionsPrefix((*ctx)->ksp,"ts_monitor_sp_eig_");CHKERRQ(ierr); /* this is wrong, used use also prefix from the TS */
  ierr = KSPSetType((*ctx)->ksp,KSPGMRES);CHKERRQ(ierr);
  ierr = KSPGMRESSetRestart((*ctx)->ksp,200);CHKERRQ(ierr);
  ierr = KSPSetTolerances((*ctx)->ksp,1.e-10,PETSC_DEFAULT,PETSC_DEFAULT,200);CHKERRQ(ierr);
  ierr = KSPSetComputeSingularValues((*ctx)->ksp,PETSC_TRUE);CHKERRQ(ierr);
  ierr = KSPSetFromOptions((*ctx)->ksp);CHKERRQ(ierr);
  ierr = KSPGetPC((*ctx)->ksp,&pc);CHKERRQ(ierr);
  ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr);

  (*ctx)->howoften          = howoften;
  (*ctx)->computeexplicitly = PETSC_FALSE;

  ierr = PetscOptionsGetBool(NULL,"-ts_monitor_sp_eig_explicitly",&(*ctx)->computeexplicitly,NULL);CHKERRQ(ierr);

  (*ctx)->comm = comm;
  (*ctx)->xmin = -2.1;
  (*ctx)->xmax = 1.1;
  (*ctx)->ymin = -1.1;
  (*ctx)->ymax = 1.1;
  PetscFunctionReturn(0);
}

  /// Do what is necessary to build this instance
  void p_build(const std::string& option_prefix)
  {
    PetscErrorCode ierr;
    try  {
      parallel::Communicator comm(this->communicator());
      if (this->p_doSerial) {
        comm = this->communicator().self();
      }
      ierr = KSPCreate(comm, &p_KSP); CHKERRXX(ierr);
      if (!this->p_guessZero) {
        ierr = KSPSetInitialGuessNonzero(p_KSP,PETSC_TRUE); CHKERRXX(ierr); 
      } else {
        ierr = KSPSetInitialGuessNonzero(p_KSP,PETSC_FALSE); CHKERRXX(ierr); 
      }
      ierr = KSPSetOptionsPrefix(p_KSP, option_prefix.c_str()); CHKERRXX(ierr);
      PC pc;
      ierr = KSPGetPC(p_KSP, &pc); CHKERRXX(ierr);
      ierr = PCSetOptionsPrefix(pc, option_prefix.c_str()); CHKERRXX(ierr);

      ierr = KSPSetTolerances(p_KSP, 
                              LinearSolverImplementation<T, I>::p_relativeTolerance, 
                              LinearSolverImplementation<T, I>::p_solutionTolerance, 
                              PETSC_DEFAULT,
                              LinearSolverImplementation<T, I>::p_maxIterations); CHKERRXX(ierr);

      ierr = KSPSetFromOptions(p_KSP);CHKERRXX(ierr);
    } catch (const PETSC_EXCEPTION_TYPE& e) {
      throw PETScException(ierr, e);
    }
  }  

	/*! \brief try to solve the system with block jacobi pre-conditioner
	 *
	 *
	 *
	 */
	void try_solve_complex_bj(Mat & A_, const Vec & b_, Vec & x_)
	{
		PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol,abstol,dtol,5));
		PETSC_SAFE_CALL(KSPSetConvergenceTest(ksp,KSPConvergedSkip,NULL,NULL));

		solve_complex(A_,b_,x_);
	}

PetscErrorCode SolvePressure(UserContext *uc)
{
  PetscErrorCode ierr;
  
  PetscFunctionBegin;
  PetscLogEventBegin(EVENT_SolvePressure,0,0,0,0);

  ierr = KSPCreate(PETSC_COMM_WORLD, &uc->ksp); CHKERRQ(ierr);
  ierr = KSPSetType(uc->ksp,KSPCG); CHKERRQ(ierr);
  ierr = KSPSetOperators(uc->ksp,uc->A,uc->A,SAME_NONZERO_PATTERN); CHKERRQ(ierr);
  ierr = KSPSetFromOptions(uc->ksp); CHKERRQ(ierr);
  ierr = KSPSetTolerances(uc->ksp, 0.000001, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT); CHKERRQ(ierr);
  
  ierr = KSPSetType( uc->ksp, KSPPREONLY); CHKERRQ(ierr);
  PC pc;
  ierr = KSPGetPC(uc->ksp, &pc); CHKERRQ(ierr);
  ierr = PCSetType(pc, PCLU); CHKERRQ(ierr);
  
  ierr = KSPSolve(uc->ksp,uc->b,uc->p);CHKERRQ(ierr);
  
  KSPConvergedReason reason;
  KSPGetConvergedReason(uc->ksp,&reason);
  PetscPrintf(PETSC_COMM_WORLD,"Pressure KSPConvergedReason: %D\n", reason);
  // if( reason < 0 )    SETERRQ(PETSC_ERR_CONV_FAILED, "Exiting: Failed to converge\n");
  
  PetscLogEventEnd(EVENT_SolvePressure,0,0,0,0);
  PetscFunctionReturn(0);
}

static PetscErrorCode KSPChebyshevEstEigSet_Chebyshev(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
{
  KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  if (a != 0.0 || b != 0.0 || c != 0.0 || d != 0.0) {
    if (!cheb->kspest) { /* should this block of code be moved to KSPSetUp_Chebyshev()? */
      ierr = KSPCreate(PetscObjectComm((PetscObject)ksp),&cheb->kspest);CHKERRQ(ierr);
      ierr = PetscObjectIncrementTabLevel((PetscObject)cheb->kspest,(PetscObject)ksp,1);CHKERRQ(ierr);
      ierr = KSPSetOptionsPrefix(cheb->kspest,((PetscObject)ksp)->prefix);CHKERRQ(ierr);
      ierr = KSPAppendOptionsPrefix(cheb->kspest,"esteig_");CHKERRQ(ierr);
      ierr = KSPSetSkipPCSetFromOptions(cheb->kspest,PETSC_TRUE);CHKERRQ(ierr);

      ierr = KSPSetPC(cheb->kspest,ksp->pc);CHKERRQ(ierr);

      ierr = KSPSetComputeEigenvalues(cheb->kspest,PETSC_TRUE);CHKERRQ(ierr);

      /* We cannot turn off convergence testing because GMRES will break down if you attempt to keep iterating after a zero norm is obtained */
      ierr = KSPSetTolerances(cheb->kspest,1.e-12,PETSC_DEFAULT,PETSC_DEFAULT,cheb->eststeps);CHKERRQ(ierr);
    }
    if (a >= 0) cheb->tform[0] = a;
    if (b >= 0) cheb->tform[1] = b;
    if (c >= 0) cheb->tform[2] = c;
    if (d >= 0) cheb->tform[3] = d;
    cheb->amatid    = 0;
    cheb->pmatid    = 0;
    cheb->amatstate = -1;
    cheb->pmatstate = -1;
  } else {
    ierr = KSPDestroy(&cheb->kspest);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

/*!
  \brief KSP and PC type

 KSPRICHARDSON "richardson"
 KSPCHEBYCHEV  "chebychev"
 KSPCG         "cg"
 KSPCGNE       "cgne"
 KSPNASH       "nash"
 KSPSTCG       "stcg"
 KSPGLTR       "gltr"
 KSPGMRES      "gmres"
 KSPFGMRES     "fgmres" 
 KSPLGMRES     "lgmres"
 KSPDGMRES     "dgmres"
 KSPTCQMR      "tcqmr"
 KSPBCGS       "bcgs"
 KSPIBCGS        "ibcgs"
 KSPBCGSL        "bcgsl"
 KSPCGS        "cgs"
 KSPTFQMR      "tfqmr"
 KSPCR         "cr"
 KSPLSQR       "lsqr"
 KSPPREONLY    "preonly"
 KSPQCG        "qcg"
 KSPBICG       "bicg"
 KSPMINRES     "minres"
 KSPSYMMLQ     "symmlq"
 KSPLCD        "lcd"
 KSPPYTHON     "python"
 KSPBROYDEN    "broyden"
 KSPGCR        "gcr"
 KSPNGMRES     "ngmres"
 KSPSPECEST    "specest"

 PCNONE            "none"
 PCJACOBI          "jacobi"
 PCSOR             "sor"
 PCLU              "lu"
 PCSHELL           "shell"
 PCBJACOBI         "bjacobi"
 PCMG              "mg"
 PCEISENSTAT       "eisenstat"
 PCILU             "ilu"
 PCICC             "icc"
 PCASM             "asm"
 PCGASM            "gasm"
 PCKSP             "ksp"
 PCCOMPOSITE       "composite"
 PCREDUNDANT       "redundant"
 PCSPAI            "spai"
 PCNN              "nn"
 PCCHOLESKY        "cholesky"
 PCPBJACOBI        "pbjacobi"
 PCMAT             "mat"
 PCHYPRE           "hypre"
 PCPARMS           "parms"
 PCFIELDSPLIT      "fieldsplit"
 PCTFS             "tfs"
 PCML              "ml"
 PCPROMETHEUS      "prometheus"
 PCGALERKIN        "galerkin"
 PCEXOTIC          "exotic"
 PCHMPI            "hmpi"
 PCSUPPORTGRAPH    "supportgraph"
 PCASA             "asa"
 PCCP              "cp"
 PCBFBT            "bfbt"
 PCLSC             "lsc"
 PCPYTHON          "python"
 PCPFMG            "pfmg"
 PCSYSPFMG         "syspfmg"
 PCREDISTRIBUTE    "redistribute"
 PCSACUSP          "sacusp"
 PCSACUSPPOLY      "sacusppoly"
 PCBICGSTABCUSP    "bicgstabcusp"
 PCSVD             "svd"
 PCAINVCUSP        "ainvcusp"
 PCGAMG            "gamg"

*/
void PETScLinearSolver::Config(const PetscReal tol, const PetscInt maxits, const KSPType lsol,
                               const PCType prec_type, const std::string &prefix)
{
   ltolerance = tol;
   sol_type = lsol;
   pc_type = prec_type; 

   KSPCreate(PETSC_COMM_WORLD,&lsolver);

#if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR > 4)
   KSPSetOperators(lsolver, A, A);
#else
   KSPSetOperators(lsolver, A, A, DIFFERENT_NONZERO_PATTERN);
#endif

   KSPSetType(lsolver,lsol);

   KSPGetPC(lsolver, &prec);
   PCSetType(prec, prec_type); //  PCJACOBI); //PCNONE);
   KSPSetTolerances(lsolver,ltolerance, PETSC_DEFAULT, PETSC_DEFAULT, maxits);
//.........这里部分代码省略.........

static PetscErrorCode MatDuplicate_LMVM(Mat B, MatDuplicateOption op, Mat *mat)
{
  Mat_LMVM          *bctx = (Mat_LMVM*)B->data;
  Mat_LMVM          *mctx;
  PetscErrorCode    ierr;
  MatType           lmvmType;
  Mat               A;

  PetscFunctionBegin;
  ierr = MatGetType(B, &lmvmType);CHKERRQ(ierr);
  ierr = MatCreate(PetscObjectComm((PetscObject)B), mat);CHKERRQ(ierr);
  ierr = MatSetType(*mat, lmvmType);CHKERRQ(ierr);

  A = *mat;
  mctx = (Mat_LMVM*)A->data;
  mctx->m = bctx->m;
  mctx->ksp_max_it = bctx->ksp_max_it;
  mctx->ksp_rtol = bctx->ksp_rtol;
  mctx->ksp_atol = bctx->ksp_atol;
  mctx->shift = bctx->shift;
  ierr = KSPSetTolerances(mctx->J0ksp, mctx->ksp_rtol, mctx->ksp_atol, PETSC_DEFAULT, mctx->ksp_max_it);CHKERRQ(ierr);

  ierr = MatLMVMAllocate(*mat, bctx->Xprev, bctx->Fprev);CHKERRQ(ierr);
  if (op == MAT_COPY_VALUES) {
    ierr = MatCopy(B, *mat, SAME_NONZERO_PATTERN);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

  void PetscPreconditioner::set_petsc_subpreconditioner_type(const PCType type, PC& pc)  {

    int ierr;
    KSP* subksps;
    int nlocal;

    ierr = PCASMGetSubKSP(pc, &nlocal, PETSC_NULL, &subksps);
    CHKERRABORT(MPI_COMM_WORLD, ierr);

    PetscReal epsilon = 1.e-16;

    for(int i = 0; i < nlocal; i++) {
      PC subpc;

      ierr = KSPGetPC(subksps[i], &subpc);
      CHKERRABORT(MPI_COMM_WORLD, ierr);

      ierr = KSPSetTolerances(subksps[i], PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, 1);
      CHKERRABORT(MPI_COMM_WORLD, ierr);

      ierr = KSPSetFromOptions(subksps[i]);
      CHKERRABORT(MPI_COMM_WORLD, ierr);

      ierr = PCSetType(subpc, type);
      CHKERRABORT(MPI_COMM_WORLD, ierr);

      ierr = PCFactorSetZeroPivot(subpc, epsilon);
      CHKERRABORT(MPI_COMM_WORLD, ierr);

      ierr = PCFactorSetShiftType(subpc, MAT_SHIFT_NONZERO);
      CHKERRABORT(MPI_COMM_WORLD, ierr);
    }
  }

/*@

 KSPFGMRESModifyPCKSP - modifies the attributes of the
     GMRES preconditioner.  It serves as an example (not as something
     useful!)

  Input Parameters:
+    ksp - the ksp context being used.
.    total_its     - the total number of FGMRES iterations that have occurred.
.    loc_its       - the number of FGMRES iterations since last restart.
.    res_norm      - the current residual norm.
-    dummy         - context, not used here

   Level: intermediate

   Contributed by Allison Baker

 This could be used as a template!

.seealso: KSPFGMRESSetModifyPC(), KSPFGMRESModifyPCKSP()

@*/
PetscErrorCode  KSPFGMRESModifyPCKSP(KSP ksp,PetscInt total_its,PetscInt loc_its,PetscReal res_norm,void *dummy)
{
  PC             pc;
  PetscErrorCode ierr;
  PetscInt       maxits;
  KSP            sub_ksp;
  PetscReal      rtol,abstol,dtol;
  PetscBool      isksp;

  PetscFunctionBegin;
  ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);

  ierr = PetscObjectTypeCompare((PetscObject)pc,PCKSP,&isksp);CHKERRQ(ierr);
  if (isksp) {
    ierr = PCKSPGetKSP(pc,&sub_ksp);CHKERRQ(ierr);

    /* note that at this point you could check the type of KSP with KSPGetType() */

    /* Now we can use functions such as KSPGMRESSetRestart() or
      KSPGMRESSetOrthogonalization() or KSPSetTolerances() */

    ierr = KSPGetTolerances(sub_ksp,&rtol,&abstol,&dtol,&maxits);CHKERRQ(ierr);
    if (!loc_its) {
      rtol = .1;
    } else {
      rtol *= .9;
    }
    ierr = KSPSetTolerances(sub_ksp,rtol,abstol,dtol,maxits);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

  void PetscDMNonlinearSolver<T>::init()
  {
    PetscErrorCode ierr;
    DM dm;
    this->PetscNonlinearSolver<T>::init();

    // Attaching a DM with the function and Jacobian callbacks to SNES.
    ierr = DMCreateLibMesh(libMesh::COMM_WORLD, this->system(), &dm); CHKERRABORT(libMesh::COMM_WORLD, ierr);
    ierr = DMSetFromOptions(dm);               CHKERRABORT(libMesh::COMM_WORLD, ierr);
    ierr = DMSetUp(dm);                        CHKERRABORT(libMesh::COMM_WORLD, ierr);
    ierr = SNESSetDM(this->_snes, dm);         CHKERRABORT(libMesh::COMM_WORLD, ierr);
    // SNES now owns the reference to dm.
    ierr = DMDestroy(&dm);                     CHKERRABORT(libMesh::COMM_WORLD, ierr);
    KSP ksp;
    ierr = SNESGetKSP (this->_snes, &ksp);     CHKERRABORT(libMesh::COMM_WORLD,ierr);

    // Set the tolerances for the iterative solver.  Use the user-supplied
    // tolerance for the relative residual & leave the others at default values
    ierr = KSPSetTolerances (ksp, this->initial_linear_tolerance, PETSC_DEFAULT,PETSC_DEFAULT, this->max_linear_iterations); CHKERRABORT(libMesh::COMM_WORLD,ierr);

    // Set the tolerances for the non-linear solver.
    ierr = SNESSetTolerances(this->_snes,
			     this->absolute_residual_tolerance,
			     this->relative_residual_tolerance,
			     this->absolute_step_tolerance,
			     this->max_nonlinear_iterations,
			     this->max_function_evaluations);
    CHKERRABORT(libMesh::COMM_WORLD,ierr);

    //Pull in command-line options
    KSPSetFromOptions(ksp);
    SNESSetFromOptions(this->_snes);
  }

/*
   PostSetSubKSP - Optional user-defined routine that reset SubKSP options when hierarchical bjacobi PC is used
   e.g,
     mpiexec -n 8 ./ex3 -nox -n 10000 -ksp_type fgmres -pc_type bjacobi -pc_bjacobi_blocks 4 -sub_ksp_type gmres -sub_ksp_max_it 3 -post_setsubksp -sub_ksp_rtol 1.e-16
   Set by SNESLineSearchSetPostCheck().

   Input Parameters:
   linesearch - the LineSearch context
   xcurrent - current solution
   y - search direction and length
   x    - the new candidate iterate

   Output Parameters:
   y    - proposed step (search direction and length) (possibly changed)
   x    - current iterate (possibly modified)

 */
PetscErrorCode PostSetSubKSP(SNESLineSearch linesearch,Vec xcurrent,Vec y,Vec x,PetscBool  *changed_y,PetscBool  *changed_x, void * ctx)
{
  PetscErrorCode ierr;
  SetSubKSPCtx   *check;
  PetscInt       iter,its,sub_its,maxit;
  KSP            ksp,sub_ksp,*sub_ksps;
  PC             pc;
  PetscReal      ksp_ratio;
  SNES           snes;

  PetscFunctionBeginUser;
  ierr    = SNESLineSearchGetSNES(linesearch, &snes);CHKERRQ(ierr);
  check   = (SetSubKSPCtx*)ctx;
  ierr    = SNESGetIterationNumber(snes,&iter);CHKERRQ(ierr);
  ierr    = SNESGetKSP(snes,&ksp);CHKERRQ(ierr);
  ierr    = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
  ierr    = PCBJacobiGetSubKSP(pc,NULL,NULL,&sub_ksps);CHKERRQ(ierr);
  sub_ksp = sub_ksps[0];
  ierr    = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);      /* outer KSP iteration number */
  ierr    = KSPGetIterationNumber(sub_ksp,&sub_its);CHKERRQ(ierr); /* inner KSP iteration number */

  if (iter) {
    ierr      = PetscPrintf(PETSC_COMM_WORLD,"    ...PostCheck snes iteration %D, ksp_it %d %d, subksp_it %d\n",iter,check->its0,its,sub_its);CHKERRQ(ierr);
    ksp_ratio = ((PetscReal)(its))/check->its0;
    maxit     = (PetscInt)(ksp_ratio*sub_its + 0.5);
    if (maxit < 2) maxit = 2;
    ierr = KSPSetTolerances(sub_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,maxit);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"    ...ksp_ratio %g, new maxit %d\n\n",ksp_ratio,maxit);CHKERRQ(ierr);
  }
  check->its0 = its; /* save current outer KSP iteration number */
  PetscFunctionReturn(0);
}

void PETSc::SetTol(double val)
{
	PetscInt maxits;
	double rtol, atol, dtol;
	
	KSPGetTolerances(ksp, &rtol, &atol, &dtol, &maxits);
	ierr = KSPSetTolerances(ksp, val, atol, dtol, maxits);
}

void MultigridSolver_SetAbsoluteTolerance( void* matrixSolver, double tolerance ) {
	MultigridSolver*	self = (MultigridSolver*)matrixSolver;

	assert( self && Stg_CheckType( self, MultigridSolver ) );

	//MatrixSolver_SetAbsoluteTolerance( self->outerSolver, tolerance );
	KSPSetTolerances( self->outerSolver->ksp, PETSC_DEFAULT, tolerance, PETSC_DEFAULT, PETSC_DEFAULT );

	self->solversChanged = True;
}

bool
IBImplicitModHelmholtzPETScLevelSolver::solveSystem(
    SAMRAIVectorReal<NDIM,double>& x,
    SAMRAIVectorReal<NDIM,double>& b)
{
    IBAMR_TIMER_START(t_solve_system);

    int ierr;

    if (d_enable_logging) plog << d_object_name << "::solveSystem():" << std::endl;

    // Initialize the solver, when necessary.
    const bool deallocate_after_solve = !d_is_initialized;
    if (deallocate_after_solve) initializeSolverState(x,b);
#if 0 // XXXX
    // Configure solver.
    ierr = KSPSetTolerances(d_petsc_ksp, d_rel_residual_tol, d_abs_residual_tol, PETSC_DEFAULT, d_max_iterations); IBTK_CHKERRQ(ierr);
    ierr = KSPSetInitialGuessNonzero(d_petsc_ksp, d_initial_guess_nonzero ? PETSC_TRUE : PETSC_FALSE); IBTK_CHKERRQ(ierr);
#endif

    // Solve the system.
    Pointer<PatchLevel<NDIM> > patch_level = d_hierarchy->getPatchLevel(d_level_num);
    const int x_idx = x.getComponentDescriptorIndex(0);
    Pointer<SideVariable<NDIM,double> > x_var = x.getComponentVariable(0);
    const int b_idx = b.getComponentDescriptorIndex(0);
    Pointer<SideVariable<NDIM,double> > b_var = b.getComponentVariable(0);
    if (d_initial_guess_nonzero) PETScVecUtilities::copyToPatchLevelVec(d_petsc_x, x_idx, x_var, patch_level);
    PETScVecUtilities::copyToPatchLevelVec(d_petsc_b, b_idx, b_var, patch_level);
    PETScVecUtilities::constrainPatchLevelVec(d_petsc_b, d_dof_index_idx, d_dof_index_var, patch_level, d_dof_index_fill);
    ierr = KSPSolve(d_petsc_ksp, d_petsc_b, d_petsc_x); IBTK_CHKERRQ(ierr);
    PETScVecUtilities::copyFromPatchLevelVec(d_petsc_x, x_idx, x_var, patch_level);

    typedef SideDataSynchronization::SynchronizationTransactionComponent SynchronizationTransactionComponent; // XXXX
    SynchronizationTransactionComponent x_synch_transaction = SynchronizationTransactionComponent(x_idx, "CONSERVATIVE_COARSEN");
    Pointer<SideDataSynchronization> side_synch_op = new SideDataSynchronization();
    side_synch_op->initializeOperatorState(x_synch_transaction, x.getPatchHierarchy());
    side_synch_op->synchronizeData(0.0);

    // Log solver info.
    KSPConvergedReason reason;
    ierr = KSPGetConvergedReason(d_petsc_ksp, &reason); IBTK_CHKERRQ(ierr);
    const bool converged = reason > 0;
    if (d_enable_logging)
    {
        plog << d_object_name << "::solveSystem(): solver " << (converged ? "converged" : "diverged") << "\n"
             << "iterations = " << d_current_its << "\n"
             << "residual norm = " << d_current_residual_norm << std::endl;
    }

    // Deallocate the solver, when necessary.
    if (deallocate_after_solve) deallocateSolverState();

    IBAMR_TIMER_STOP(t_solve_system);
    return converged;
}// solveSystem

/* This function returns the suitable solver for the velocity linear system */
KSP Solver_get_velocity_solver(Mat lhs, Parameters *params) {

    KSP solver;
    PC pc;
    double rtol;
    int ierr;

    rtol = params->resmax;

    ierr = KSPCreate(PETSC_COMM_WORLD, &solver); PETScErrAct(ierr);
    ierr = KSPSetOperators(solver, lhs, lhs, SAME_NONZERO_PATTERN); PETScErrAct(ierr);

    //ierr = KSPSetType(solver, KSPGMRES); PETScErrAct(ierr);
    ierr = KSPSetType(solver, KSPIBCGS); PETScErrAct(ierr);

    ierr = KSPSetTolerances(solver, rtol, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT); PETScErrAct(ierr);

    ierr = KSPGetPC(solver, &pc); PETScErrAct(ierr);
    ierr = PCSetType(pc, PCNONE); PETScErrAct(ierr);
    //ierr = PCSetType(pc, PCJACOBI); PETScErrAct(ierr);
    //ierr = PCSetType(pc, PCNONE); PETScErrAct(ierr);

    //ierr = PCSetType(pc, PCASM); PETScErrAct(ierr);

    //ierr = PCSetType(pc, PCBJACOBI); PETScErrAct(ierr);

    //ierr = PCSetType(pc, PCBJACOBI); PETScErrAct(ierr);
    //ierr = PCSetType(pc, PCLU); PETScErrAct(ierr);

    //ierr = PCSetType(pc, PCEISENSTAT); PETScErrAct(ierr);
    //ierr = PCSetType(pc, PCSOR); PETScErrAct(ierr);
    //ierr = PCSetType(pc, PCJACOBI); PETScErrAct(ierr);
    //ierr = PCSetType(pc, PCNONE); PETScErrAct(ierr);

    /*
 ierr = KSPGetPC(solver, &pc); PETScErrAct(ierr);
 ierr = PCSetType(pc, PCILU); PETScErrAct(ierr);
 ierr = PCFactorSetLevels(pc, 3); PETScErrAct(ierr);

 //ierr = PCFactorSetUseDropTolerance(pc, 1e-3, .1, 50); PETScErrAct(ierr);
 //ierr = PCFactorSetFill(pc, 30.7648); PETScErrAct(ierr);

 ierr = PCFactorSetReuseOrdering(pc, PETSC_TRUE); PETScErrAct(ierr);
 ierr = PCFactorSetReuseFill(pc, PETSC_TRUE); PETScErrAct(ierr);
 ierr = PCFactorSetAllowDiagonalFill(pc); PETScErrAct(ierr);

 //ierr = PCFactorSetUseInPlace(pc); PETScErrAct(ierr);
*/
    ierr = KSPSetInitialGuessNonzero(solver, PETSC_TRUE); PETScErrAct(ierr);
    ierr = KSPSetFromOptions(solver); PETScErrAct(ierr);


    return solver;
}

void createOuterKsp(OuterContext* ctx) {
  KSPCreate((ctx->data)->commAll, &(ctx->outerKsp));
  PetscObjectIncrementTabLevel((PetscObject)(ctx->outerKsp), PETSC_NULL, 0);
  KSPSetType(ctx->outerKsp, KSPFGMRES);
  KSPSetTolerances(ctx->outerKsp, 1.0e-12, 1.0e-12, PETSC_DEFAULT, 50);
  KSPSetPC(ctx->outerKsp, ctx->outerPC);
  KSPSetOptionsPrefix(ctx->outerKsp, "outer_");
  KSPSetFromOptions(ctx->outerKsp);
  KSPSetOperators(ctx->outerKsp, ctx->outerMat,
      ctx->outerMat, SAME_NONZERO_PATTERN);
  KSPSetUp(ctx->outerKsp);
}