C++ KSPSetFromOptions函数代码示例

void PETSc::Solve_withPureNeumann_GMRES(void)
{
    	ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
  	ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
  	
  	ierr = VecAssemblyBegin(x);
  	ierr = VecAssemblyEnd(x);
  	
  	ierr = VecAssemblyBegin(b);
  	ierr = VecAssemblyEnd(b);
  	
	
	MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,PETSC_NULL,&nullsp);
        KSPSetNullSpace(ksp,nullsp);
	MatNullSpaceRemove(nullsp,b,PETSC_NULL);

	
        KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);
        
	KSPSetType(ksp,KSPGMRES);

	//KSPGetPC(ksp, &pc);
	//PCSetType(pc, PCASM);
        KSPSetFromOptions(ksp);
        KSPSetUp(ksp);
	//start_clock("Before Petsc Solve in pure neumann solver");
        KSPSolve(ksp,b,x);
	//stop_clock("After Petsc Solve in pure neumann solver");
}

int main(int argc,char **argv)
{
  KSP            ksp;
  DM             da;
  UserContext    user;
  PetscInt       bc;
  PetscErrorCode ierr;

  PetscInitialize(&argc,&argv,(char *)0,help);
  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
  ierr = DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-11,-11,PETSC_DECIDE,PETSC_DECIDE,1,1,PETSC_NULL,PETSC_NULL,&da);CHKERRQ(ierr);
  ierr = KSPSetDM(ksp,(DM)da);
  ierr = DMSetApplicationContext(da,&user);CHKERRQ(ierr);

  user.uu = 1.0;
  user.tt = 1.0;
  bc   = (PetscInt)NEUMANN; // Use Neumann Boundary Conditions
  user.bcType = (BCType)bc;


  ierr = KSPSetComputeRHS(ksp,ComputeRHS,&user);CHKERRQ(ierr);
  ierr = KSPSetComputeOperators(ksp,ComputeJacobian,&user);CHKERRQ(ierr);
  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
  ierr = KSPSolve(ksp,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);

  ierr = DMDestroy(&da);CHKERRQ(ierr);
  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
  ierr = PetscFinalize();CHKERRQ(ierr);
  return 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);
    }
  }  

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);
}

int main(int argc, char **argv)
{
    Mat            A;
    KSP            ksp;
    DM             shell;
    Vec            *left, *right;
    MPI_Comm       c;
    PetscErrorCode ierr;

    ierr = PetscInitialize(&argc, &argv, NULL, NULL);if (ierr) return ierr;
    c = PETSC_COMM_WORLD;

    ierr = MatCreate(c, &A); CHKERRQ(ierr);
    ierr = MatSetSizes(A, 1, 1, PETSC_DECIDE, PETSC_DECIDE); CHKERRQ(ierr);
    ierr = MatSetFromOptions(A); CHKERRQ(ierr);
    ierr = MatSetUp(A); CHKERRQ(ierr);
    ierr = KSPCreate(c, &ksp); CHKERRQ(ierr);
    ierr = KSPSetOperators(ksp, A, A); CHKERRQ(ierr);
    ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr);
    ierr = DMShellCreate(c, &shell); CHKERRQ(ierr);
    ierr = DMSetFromOptions(shell); CHKERRQ(ierr);
    ierr = DMSetUp(shell); CHKERRQ(ierr);
    ierr = KSPSetDM(ksp, shell); CHKERRQ(ierr);

    ierr = KSPCreateVecs(ksp, 1, &right, 1, &left); CHKERRQ(ierr);
    ierr = VecView(right[0], PETSC_VIEWER_STDOUT_(c));CHKERRQ(ierr);
    ierr = VecDestroyVecs(1,&right); CHKERRQ(ierr);
    ierr = VecDestroyVecs(1,&left); CHKERRQ(ierr);

    ierr = DMDestroy(&shell); CHKERRQ(ierr);
    ierr = KSPDestroy(&ksp); CHKERRQ(ierr);
    ierr = MatDestroy(&A); CHKERRQ(ierr);
    PetscFinalize();
    return 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);
    }
  }

int main(int argc,char **args)
{
  KSP            ksp;             /* linear solver context */
  Mat            A;            /* matrix */
  Vec            x,b;          /* approx solution, RHS, exact solution */
  PetscViewer    fd;               /* viewer */
  PetscErrorCode ierr;

  PetscInitialize(&argc,&args,(char *)0,help);
  fd = PETSC_VIEWER_SOCKET_WORLD;

  ierr = VecCreate(PETSC_COMM_WORLD,&b);CHKERRQ(ierr);
  ierr = VecLoad(b,fd);CHKERRQ(ierr);
  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = MatLoad(A,fd);CHKERRQ(ierr);
  ierr = VecDuplicate(b,&x);CHKERRQ(ierr);

  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
  ierr = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
  ierr = KSPSetUp(ksp);CHKERRQ(ierr);
  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
  ierr = VecView(x,fd);CHKERRQ(ierr);
  ierr = MatDestroy(&A);CHKERRQ(ierr);
  ierr = VecDestroy(&b);CHKERRQ(ierr);
  ierr = VecDestroy(&x);CHKERRQ(ierr);
  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);

  ierr = PetscFinalize();
  return 0;
}

int main(int argc,char **argv)
{
  PetscErrorCode ierr;
  KSP            ksp;
  DM             da,shell;
  PetscInt       levels;

  PetscInitialize(&argc,&argv,(char*)0,help);

  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
  ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,-129,1,1,0,&da);CHKERRQ(ierr);
  ierr = MyDMShellCreate(PETSC_COMM_WORLD,da,&shell);CHKERRQ(ierr);
  /* these two lines are not needed but allow PCMG to automatically know how many multigrid levels the user wants */
  ierr = DMGetRefineLevel(da,&levels);CHKERRQ(ierr);
  ierr = DMSetRefineLevel(shell,levels);CHKERRQ(ierr);

  ierr = KSPSetDM(ksp,shell);CHKERRQ(ierr);
  ierr = KSPSetComputeRHS(ksp,ComputeRHS,NULL);CHKERRQ(ierr);
  ierr = KSPSetComputeOperators(ksp,ComputeMatrix,NULL);CHKERRQ(ierr);
  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
  ierr = KSPSolve(ksp,NULL,NULL);CHKERRQ(ierr);

  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
  ierr = DMDestroy(&shell);CHKERRQ(ierr);
  ierr = DMDestroy(&da);CHKERRQ(ierr);
  ierr = PetscFinalize();

  return 0;
}

static PetscErrorCode TaoSetFromOptions_BQNLS(PetscOptionItems *PetscOptionsObject,Tao tao)
{
  TAO_BNK        *bnk = (TAO_BNK *)tao->data;
  TAO_BQNK       *bqnk = (TAO_BQNK*)bnk->ctx;
  PetscErrorCode ierr;
  KSPType        ksp_type;
  PetscBool      is_spd;

  PetscFunctionBegin;
  ierr = PetscOptionsHead(PetscOptionsObject,"Quasi-Newton-Krylov method for bound constrained optimization");CHKERRQ(ierr);
  ierr = PetscOptionsEList("-tao_bqnls_as_type", "active set estimation method", "", BNK_AS, BNK_AS_TYPES, BNK_AS[bnk->as_type], &bnk->as_type, 0);CHKERRQ(ierr);
  ierr = PetscOptionsReal("-tao_bqnls_epsilon", "(developer) tolerance used when computing actual and predicted reduction", "", bnk->epsilon, &bnk->epsilon,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal("-tao_bqnls_as_tol", "(developer) initial tolerance used when estimating actively bounded variables", "", bnk->as_tol, &bnk->as_tol,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal("-tao_bqnls_as_step", "(developer) step length used when estimating actively bounded variables", "", bnk->as_step, &bnk->as_step,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsInt("-tao_bqnls_max_cg_its", "number of BNCG iterations to take for each Newton step", "", bnk->max_cg_its, &bnk->max_cg_its,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsTail();CHKERRQ(ierr);
  ierr = TaoSetFromOptions(bnk->bncg);CHKERRQ(ierr);
  ierr = TaoLineSearchSetFromOptions(tao->linesearch);CHKERRQ(ierr);
  ierr = KSPSetFromOptions(tao->ksp);CHKERRQ(ierr);
  ierr = KSPGetType(tao->ksp,&ksp_type);CHKERRQ(ierr);
  bnk->is_nash = bnk->is_gltr = bnk->is_stcg = PETSC_FALSE;
  ierr = MatSetFromOptions(bqnk->B);CHKERRQ(ierr);
  ierr = MatGetOption(bqnk->B, MAT_SPD, &is_spd);CHKERRQ(ierr);
  if (!is_spd) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix must be symmetric positive-definite");
  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);
}

//.........这里部分代码省略.........
 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);

   if( !prefix.empty() )
   {
       KSPSetOptionsPrefix(lsolver, prefix.c_str());
       PCSetOptionsPrefix(prec, prefix.c_str());
   }

   KSPSetFromOptions(lsolver);
}

static PetscErrorCode PCSetUp_KSP(PC pc)
{
  PetscErrorCode ierr;
  PC_KSP         *jac = (PC_KSP*)pc->data;
  Mat            mat;
  PetscBool      A;

  PetscFunctionBegin;
  if (!jac->ksp) {ierr = PCKSPCreateKSP_KSP(pc);CHKERRQ(ierr);}
  ierr = KSPSetFromOptions(jac->ksp);CHKERRQ(ierr);
  if (jac->use_true_matrix) mat = pc->mat;
  else                      mat = pc->pmat;

  ierr = KSPGetOperatorsSet(jac->ksp,&A,PETSC_NULL);CHKERRQ(ierr);
  if (!A) {
    ierr = KSPSetOperators(jac->ksp,mat,pc->pmat,pc->flag);CHKERRQ(ierr);
  } else if (pc->flag != SAME_PRECONDITIONER) {
    Mat Amat,Bmat;
    ierr = KSPGetOperators(jac->ksp,&Amat,&Bmat,PETSC_NULL);CHKERRQ(ierr);
    if (Amat == mat && Bmat == pc->pmat) {
      /* The user has not replaced the matrices so we are expected to forward the update. This incorrectly diagnoses
       * changed matrices at the top level as the user manually changing the inner matrices, but we have no way to
       * identify that in this context. The longer term solution is to track matrix state internally.
       */
      ierr = KSPSetOperators(jac->ksp,mat,pc->pmat,pc->flag);CHKERRQ(ierr);
    }
  }
  ierr = KSPSetUp(jac->ksp);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

int main(int argc,char **argv)
{
  PetscErrorCode ierr;
  PetscInt       i;
  KSP            ksp;
  DM             da;
  Vec            x;

  PetscInitialize(&argc,&argv,(char*)0,help);

  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
  ierr = DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC,-3,2,1,0,&da);CHKERRQ(ierr);
  ierr = KSPSetDM(ksp,da);CHKERRQ(ierr);
  ierr = KSPSetComputeRHS(ksp,ComputeRHS,NULL);CHKERRQ(ierr);
  ierr = KSPSetComputeOperators(ksp,ComputeMatrix,NULL);CHKERRQ(ierr);

  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
  ierr = DMCreateGlobalVector(da,&x);CHKERRQ(ierr);
  ierr = ComputeInitialSolution(da,x);CHKERRQ(ierr);
  ierr = DMSetApplicationContext(da,x);CHKERRQ(ierr);
  ierr = KSPSetUp(ksp);CHKERRQ(ierr);
  ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);
  for (i=0; i<10; i++) {
    ierr = KSPSolve(ksp,NULL,x);CHKERRQ(ierr);
    ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);
  }
  ierr = VecDestroy(&x);CHKERRQ(ierr);
  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
  ierr = DMDestroy(&da);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return 0;
}

static PetscErrorCode PCSetUp_LSC(PC pc)
{
  PC_LSC         *lsc = (PC_LSC*)pc->data;
  Mat            L,Lp,B,C;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = PCLSCAllocate_Private(pc);CHKERRQ(ierr);
  ierr = PetscObjectQuery((PetscObject)pc->mat,"LSC_L",(PetscObject*)&L);CHKERRQ(ierr);
  if (!L) {ierr = PetscObjectQuery((PetscObject)pc->pmat,"LSC_L",(PetscObject*)&L);CHKERRQ(ierr);}
  ierr = PetscObjectQuery((PetscObject)pc->pmat,"LSC_Lp",(PetscObject*)&Lp);CHKERRQ(ierr);
  if (!Lp) {ierr = PetscObjectQuery((PetscObject)pc->mat,"LSC_Lp",(PetscObject*)&Lp);CHKERRQ(ierr);}
  if (!L) {
    ierr = MatSchurComplementGetSubMatrices(pc->mat,NULL,NULL,&B,&C,NULL);CHKERRQ(ierr);
    if (!lsc->L) {
      ierr = MatMatMult(C,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&lsc->L);CHKERRQ(ierr);
    } else {
      ierr = MatMatMult(C,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&lsc->L);CHKERRQ(ierr);
    }
    Lp = L = lsc->L;
  }
  if (lsc->scale) {
    Mat Ap;
    ierr = MatSchurComplementGetSubMatrices(pc->mat,NULL,&Ap,NULL,NULL,NULL);CHKERRQ(ierr);
    ierr = MatGetDiagonal(Ap,lsc->scale);CHKERRQ(ierr); /* Should be the mass matrix, but we don't have plumbing for that yet */
    ierr = VecReciprocal(lsc->scale);CHKERRQ(ierr);
  }
  ierr = KSPSetOperators(lsc->kspL,L,Lp);CHKERRQ(ierr);
  ierr = KSPSetFromOptions(lsc->kspL);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);
  }