C++ KSPSetUp函数代码示例

void PETSc::Solve(void)
{
        
  //start_clock("Before Assemble matrix and vector");
	ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
  	ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
  	
  	ierr = VecAssemblyBegin(x);
  	ierr = VecAssemblyEnd(x);
  	
  	ierr = VecAssemblyBegin(b);
  	ierr = VecAssemblyEnd(b);
	//stop_clock("After Assembly matrix and vector");


        KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);
        KSPSetType(ksp,KSPBCGSL);
	KSPBCGSLSetEll(ksp,2);

	//KSPGetPC(ksp, &pc);
	//PCSetType(pc, PCJACOBI);

        KSPSetFromOptions(ksp);
        KSPSetUp(ksp);

	//start_clock("Before KSPSolve");
        KSPSolve(ksp,b,x);
	//stop_clock("After KSPSolve");
}

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

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

double SFieldSolveFor(SField sfv, double *Y, unsigned int yCount) {
    mySField sf = static_cast<mySField>(sfv);
    assert(yCount <= sf->maxN);
    assert(Y);
    assert(sf->running);
    sf->Y = Y;
    sf->curN = yCount;

    // -------------- SOLVE
    PetscErrorCode ierr;
    PetscLogDouble tic,toc;
    PetscTime(&tic);
    int pt[sf->d];

    ierr = MatZeroEntries(sf->J); CHKERRQ(ierr);
    JacobianOnD(sf->J, sf->F, 0, pt, sf);

    ierr = MatAssemblyBegin(sf->J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    ierr = MatAssemblyEnd(sf->J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);

    PetscTime(&toc);
    sf->timeAssembly += toc-tic;
    PetscTime(&tic);
    ierr = VecZeroEntries(sf->U);  CHKERRQ(ierr);
    ierr = KSPSetOperators(sf->ksp, sf->J, sf->J); CHKERRQ(ierr);
    ierr = KSPSetUp(sf->ksp);  CHKERRQ(ierr);
    ierr = KSPSolve(sf->ksp,sf->F,sf->U); CHKERRQ(ierr);
    PetscTime(&toc);
    sf->timeSolver  += toc-tic;
    return Integrate(sf->U,pt,0,sf);
}

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

/*@
  MatSchurComplementComputeExplicitOperator - Compute the Schur complement matrix explicitly

  Collective on Mat

  Input Parameter:
. M - the matrix obtained with MatCreateSchurComplement()

  Output Parameter:
. S - the Schur complement matrix

  Note: This can be expensive, so it is mainly for testing

  Level: advanced

.seealso: MatCreateSchurComplement(), MatSchurComplementUpdate()
@*/
PetscErrorCode MatSchurComplementComputeExplicitOperator(Mat M, Mat *S)
{
  Mat            B, C, D;
  KSP            ksp;
  PC             pc;
  PetscBool      isLU, isILU;
  PetscReal      fill = 2.0;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = MatSchurComplementGetSubMatrices(M, NULL, NULL, &B, &C, &D);CHKERRQ(ierr);
  ierr = MatSchurComplementGetKSP(M, &ksp);CHKERRQ(ierr);
  ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);
  ierr = PetscObjectTypeCompare((PetscObject) pc, PCLU, &isLU);CHKERRQ(ierr);
  ierr = PetscObjectTypeCompare((PetscObject) pc, PCILU, &isILU);CHKERRQ(ierr);
  if (isLU || isILU) {
    Mat       fact, Bd, AinvB, AinvBd;
    PetscReal eps = 1.0e-10;

    /* This can be sped up for banded LU */
    ierr = KSPSetUp(ksp);CHKERRQ(ierr);
    ierr = PCFactorGetMatrix(pc, &fact);CHKERRQ(ierr);
    ierr = MatConvert(B, MATDENSE, MAT_INITIAL_MATRIX, &Bd);CHKERRQ(ierr);
    ierr = MatDuplicate(Bd, MAT_DO_NOT_COPY_VALUES, &AinvBd);CHKERRQ(ierr);
    ierr = MatMatSolve(fact, Bd, AinvBd);CHKERRQ(ierr);
    ierr = MatDestroy(&Bd);CHKERRQ(ierr);
    ierr = MatChop(AinvBd, eps);CHKERRQ(ierr);
    ierr = MatConvert(AinvBd, MATAIJ, MAT_INITIAL_MATRIX, &AinvB);CHKERRQ(ierr);
    ierr = MatDestroy(&AinvBd);CHKERRQ(ierr);
    ierr = MatMatMult(C, AinvB, MAT_INITIAL_MATRIX, fill, S);CHKERRQ(ierr);
    ierr = MatDestroy(&AinvB);CHKERRQ(ierr);
  } else {
    Mat Ainvd, Ainv;

    ierr = PCComputeExplicitOperator(pc, &Ainvd);CHKERRQ(ierr);
    ierr = MatConvert(Ainvd, MATAIJ, MAT_INITIAL_MATRIX, &Ainv);CHKERRQ(ierr);
    ierr = MatDestroy(&Ainvd);CHKERRQ(ierr);
#if 0
    /* Symmetric version */
    ierr = MatPtAP(Ainv, B, MAT_INITIAL_MATRIX, fill, S);CHKERRQ(ierr);
#else
    /* Nonsymmetric version */
    ierr = MatMatMatMult(C, Ainv, B, MAT_INITIAL_MATRIX, fill, S);CHKERRQ(ierr);
#endif
    ierr = MatDestroy(&Ainv);CHKERRQ(ierr);
  }

  ierr = PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr);
  ierr = MatView(*S, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  ierr = PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

  if (D) {
    MatInfo info;

    ierr = MatGetInfo(D, MAT_GLOBAL_SUM, &info);CHKERRQ(ierr);
    if (info.nz_used) SETERRQ(PetscObjectComm((PetscObject) M), PETSC_ERR_SUP, "Not yet implemented");
  }
  PetscFunctionReturn(0);
}

static PetscErrorCode KSPComputeShifts_GMRES(KSP ksp)
{
  PetscErrorCode ierr;
  KSP_AGMRES     *agmres = (KSP_AGMRES*)(ksp->data);
  KSP            kspgmres;
  Mat            Amat, Pmat;
  MatStructure   flag;
  PetscInt       max_k = agmres->max_k;
  PC             pc;
  PetscInt       m;
  PetscScalar    *Rshift, *Ishift;


  PetscFunctionBegin;
  /* Perform one cycle of classical GMRES (with the Arnoldi process) to get the Hessenberg matrix
   We assume here that the ksp is AGMRES and that the operators for the
   linear system have been set in this ksp */
  ierr = KSPCreate(PetscObjectComm((PetscObject)ksp), &kspgmres);CHKERRQ(ierr);
  if (!ksp->pc) { ierr = KSPGetPC(ksp,&ksp->pc);CHKERRQ(ierr); }
  ierr = PCGetOperators(ksp->pc, &Amat, &Pmat);CHKERRQ(ierr);
  ierr = KSPSetOperators(kspgmres, Amat, Pmat);CHKERRQ(ierr);
  ierr = KSPSetFromOptions(kspgmres);CHKERRQ(ierr);
  ierr = PetscOptionsHasName(NULL, "-ksp_view", &flg);CHKERRQ(ierr);
  if (flag) { ierr = PetscOptionsClearValue("-ksp_view");CHKERRQ(ierr); }
  ierr = KSPSetType(kspgmres, KSPGMRES);CHKERRQ(ierr);
  ierr = KSPGMRESSetRestart(kspgmres, max_k);CHKERRQ(ierr);
  ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);
  ierr = KSPSetPC(kspgmres, pc);CHKERRQ(ierr);
  /* Copy common options */
  kspgmres->pc_side = ksp->pc_side;
  /* Setup KSP context */
  ierr = KSPSetComputeEigenvalues(kspgmres, PETSC_TRUE);CHKERRQ(ierr);
  ierr = KSPSetUp(kspgmres);CHKERRQ(ierr);

  kspgmres->max_it = max_k; /* Restrict the maximum number of iterations to one cycle of GMRES */
  kspgmres->rtol   = ksp->rtol;

  ierr = KSPSolve(kspgmres, ksp->vec_rhs, ksp->vec_sol);CHKERRQ(ierr);

  ksp->guess_zero = PETSC_FALSE;
  ksp->rnorm      = kspgmres->rnorm;
  ksp->its        = kspgmres->its;
  if (kspgmres->reason == KSP_CONVERGED_RTOL) {
    ksp->reason = KSP_CONVERGED_RTOL;
    PetscFunctionReturn(0);
  } else ksp->reason = KSP_CONVERGED_ITERATING;
  /* Now, compute the Shifts values */
  ierr = PetscMalloc2(max_k,&Rshift,max_k,&Ishift);CHKERRQ(ierr);
  ierr = KSPComputeEigenvalues(kspgmres, max_k, Rshift, Ishift, &m);CHKERRQ(ierr);
  if (m < max_k) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB, "Unable to compute the Shifts for the Newton basis");
  else {
    ierr = KSPAGMRESLejaOrdering(Rshift, Ishift, agmres->Rshift, agmres->Ishift, max_k);CHKERRQ(ierr);

    agmres->HasShifts = PETSC_TRUE;
  }
  /* Restore KSP view options */
  if (flg) { ierr = PetscOptionsSetValue("-ksp_view", "");CHKERRQ(ierr); }
  PetscFunctionReturn(0);
}

int main(int argc,char **argv)
{
    KSP            ksp;
    DM             da;
    UserContext    user;
    const char     *bcTypes[2] = {"dirichlet","neumann"};
    PetscErrorCode ierr;
    PetscInt       bc;
    Vec            b,x;

    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,-3,-3,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);
    CHKERRQ(ierr);
    ierr = DMDASetUniformCoordinates(da,0,1,0,1,0,0);
    CHKERRQ(ierr);
    ierr = DMDASetFieldName(da,0,"Pressure");
    CHKERRQ(ierr);

    ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for the inhomogeneous Poisson equation", "DMqq");
    user.rho    = 1.0;
    ierr        = PetscOptionsReal("-rho", "The conductivity", "ex29.c", user.rho, &user.rho, PETSC_NULL);
    CHKERRQ(ierr);
    user.nu     = 0.1;
    ierr        = PetscOptionsReal("-nu", "The width of the Gaussian source", "ex29.c", user.nu, &user.nu, PETSC_NULL);
    CHKERRQ(ierr);
    bc          = (PetscInt)DIRICHLET;
    ierr        = PetscOptionsEList("-bc_type","Type of boundary condition","ex29.c",bcTypes,2,bcTypes[0],&bc,PETSC_NULL);
    CHKERRQ(ierr);
    user.bcType = (BCType)bc;
    ierr = PetscOptionsEnd();

    ierr = KSPSetComputeRHS(ksp,ComputeRHS,&user);
    CHKERRQ(ierr);
    ierr = KSPSetComputeOperators(ksp,ComputeMatrix,&user);
    CHKERRQ(ierr);
    ierr = KSPSetDM(ksp,da);
    CHKERRQ(ierr);
    ierr = KSPSetFromOptions(ksp);
    CHKERRQ(ierr);
    ierr = KSPSetUp(ksp);
    CHKERRQ(ierr);
    ierr = KSPSolve(ksp,PETSC_NULL,PETSC_NULL);
    CHKERRQ(ierr);
    ierr = KSPGetSolution(ksp,&x);
    CHKERRQ(ierr);
    ierr = KSPGetRhs(ksp,&b);
    CHKERRQ(ierr);

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

    return 0;
}

static PetscErrorCode PCBDDCScalingSetUp_Deluxe_Private(PC pc)
{
  PC_BDDC                *pcbddc=(PC_BDDC*)pc->data;
  PCBDDCDeluxeScaling    deluxe_ctx=pcbddc->deluxe_ctx;
  PCBDDCSubSchurs        sub_schurs = pcbddc->sub_schurs;
  PetscErrorCode         ierr;

  PetscFunctionBegin;
  if (!sub_schurs->n_subs) {
    PetscFunctionReturn(0);
  }

  /* Create work vectors for sequential part of deluxe */
  ierr = MatCreateVecs(sub_schurs->S_Ej_all,&deluxe_ctx->seq_work1,&deluxe_ctx->seq_work2);CHKERRQ(ierr);

  /* Compute deluxe sequential scatter */
  if (sub_schurs->reuse_mumps && !sub_schurs->is_dir) {
    PCBDDCReuseMumps reuse_mumps = sub_schurs->reuse_mumps;
    ierr = PetscObjectReference((PetscObject)reuse_mumps->correction_scatter_B);CHKERRQ(ierr);
    deluxe_ctx->seq_scctx = reuse_mumps->correction_scatter_B;
  } else {
    ierr = VecScatterCreate(pcbddc->work_scaling,sub_schurs->is_Ej_all,deluxe_ctx->seq_work1,NULL,&deluxe_ctx->seq_scctx);CHKERRQ(ierr);
  }

  /* Create Mat object for deluxe scaling */
  ierr = PetscObjectReference((PetscObject)sub_schurs->S_Ej_all);CHKERRQ(ierr);
  deluxe_ctx->seq_mat = sub_schurs->S_Ej_all;
  if (sub_schurs->sum_S_Ej_all) { /* if this matrix is present, then we need to create the KSP object to invert it */
    PC               pc_temp;
    MatSolverPackage solver=NULL;
    char             ksp_prefix[256];
    size_t           len;

    ierr = KSPCreate(PETSC_COMM_SELF,&deluxe_ctx->seq_ksp);CHKERRQ(ierr);
    ierr = KSPSetOperators(deluxe_ctx->seq_ksp,sub_schurs->sum_S_Ej_all,sub_schurs->sum_S_Ej_all);CHKERRQ(ierr);
    ierr = KSPSetType(deluxe_ctx->seq_ksp,KSPPREONLY);CHKERRQ(ierr);
    ierr = KSPGetPC(deluxe_ctx->seq_ksp,&pc_temp);CHKERRQ(ierr);
    ierr = PCSetType(pc_temp,PCLU);CHKERRQ(ierr);
    ierr = KSPGetPC(pcbddc->ksp_D,&pc_temp);CHKERRQ(ierr);
    ierr = PCFactorGetMatSolverPackage(pc_temp,(const MatSolverPackage*)&solver);CHKERRQ(ierr);
    if (solver) {
      PC     new_pc;
      PCType type;

      ierr = PCGetType(pc_temp,&type);CHKERRQ(ierr);
      ierr = KSPGetPC(deluxe_ctx->seq_ksp,&new_pc);CHKERRQ(ierr);
      ierr = PCSetType(new_pc,type);CHKERRQ(ierr);
      ierr = PCFactorSetMatSolverPackage(new_pc,solver);CHKERRQ(ierr);
    }
    ierr = PetscStrlen(((PetscObject)(pcbddc->ksp_D))->prefix,&len);CHKERRQ(ierr);
    len -= 10; /* remove "dirichlet_" */
    ierr = PetscStrncpy(ksp_prefix,((PetscObject)(pcbddc->ksp_D))->prefix,len+1);CHKERRQ(ierr);
    ierr = PetscStrcat(ksp_prefix,"deluxe_");CHKERRQ(ierr);
    ierr = KSPSetOptionsPrefix(deluxe_ctx->seq_ksp,ksp_prefix);CHKERRQ(ierr);
    ierr = KSPSetFromOptions(deluxe_ctx->seq_ksp);CHKERRQ(ierr);
    ierr = KSPSetUp(deluxe_ctx->seq_ksp);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

static PetscErrorCode TaoSetup_BLMVM(Tao tao)
{
  TAO_BLMVM      *blmP = (TAO_BLMVM *)tao->data;
  PetscInt       n,N;
  PetscErrorCode ierr;
  KSP            H0ksp;

  PetscFunctionBegin;
  /* Existence of tao->solution checked in TaoSetup() */
  ierr = VecDuplicate(tao->solution,&blmP->Xold);CHKERRQ(ierr);
  ierr = VecDuplicate(tao->solution,&blmP->Gold);CHKERRQ(ierr);
  ierr = VecDuplicate(tao->solution, &blmP->unprojected_gradient);CHKERRQ(ierr);

  if (!tao->stepdirection) {
    ierr = VecDuplicate(tao->solution, &tao->stepdirection);CHKERRQ(ierr);
  }
  if (!tao->gradient) {
    ierr = VecDuplicate(tao->solution,&tao->gradient);CHKERRQ(ierr);
  }
  if (!tao->XL) {
    ierr = VecDuplicate(tao->solution,&tao->XL);CHKERRQ(ierr);
    ierr = VecSet(tao->XL,PETSC_NINFINITY);CHKERRQ(ierr);
  }
  if (!tao->XU) {
    ierr = VecDuplicate(tao->solution,&tao->XU);CHKERRQ(ierr);
    ierr = VecSet(tao->XU,PETSC_INFINITY);CHKERRQ(ierr);
  }
  /* Create matrix for the limited memory approximation */
  ierr = VecGetLocalSize(tao->solution,&n);CHKERRQ(ierr);
  ierr = VecGetSize(tao->solution,&N);CHKERRQ(ierr);
  ierr = MatCreateLMVM(((PetscObject)tao)->comm,n,N,&blmP->M);CHKERRQ(ierr);
  ierr = MatLMVMAllocateVectors(blmP->M,tao->solution);CHKERRQ(ierr);

  /* If the user has set a matrix to solve as the initial H0, set the options prefix here, and set up the KSP */
  if (blmP->H0) {
    const char *prefix;
    PC H0pc;

    ierr = MatLMVMSetH0(blmP->M, blmP->H0);CHKERRQ(ierr);
    ierr = MatLMVMGetH0KSP(blmP->M, &H0ksp);CHKERRQ(ierr);

    ierr = TaoGetOptionsPrefix(tao, &prefix);CHKERRQ(ierr);
    ierr = KSPSetOptionsPrefix(H0ksp, prefix);CHKERRQ(ierr);
    ierr = PetscObjectAppendOptionsPrefix((PetscObject)H0ksp, "tao_h0_");CHKERRQ(ierr);
    ierr = KSPGetPC(H0ksp, &H0pc);CHKERRQ(ierr);
    ierr = PetscObjectAppendOptionsPrefix((PetscObject)H0pc,  "tao_h0_");CHKERRQ(ierr);

    ierr = KSPSetFromOptions(H0ksp);CHKERRQ(ierr);
    ierr = KSPSetUp(H0ksp);CHKERRQ(ierr);
  }

  PetscFunctionReturn(0);
}

void SingleLongPipe::initialize() {
  MatCreate(PETSC_COMM_WORLD, &M);
  MatSetSizes(M, local_n, local_n, global_n, global_n);
  MatSetFromOptions(M);
  MatMPIAIJSetPreallocation(M, 2, nullptr, 2, nullptr);
  MatSeqAIJSetPreallocation(M, 2, nullptr);
  MatSetUp(M);
  KSPCreate(PETSC_COMM_WORLD, &ksp);
  KSPSetOperators(ksp, M, M);
  KSPSetFromOptions(ksp);
  KSPSetUp(ksp);
}

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

int main(int argc,char **args)
{
  KSP            ksp;
  Mat            A;
  Vec            x,b;
  PetscViewer    fd;
  char           file[PETSC_MAX_PATH_LEN];
  PetscErrorCode ierr;
  PetscBool      flg,preload = PETSC_TRUE;

  PetscInitialize(&argc,&args,(char*)0,help);
  ierr = PetscLogDefaultBegin();CHKERRQ(ierr);
  ierr = PetscOptionsGetString(NULL,NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);

  PetscPreLoadBegin(preload,"Load system");

  /*
     Load the matrix and vector; then destroy the viewer.
  */
  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = MatSetFromOptions(A);CHKERRQ(ierr);
  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);CHKERRQ(ierr);
  ierr = MatLoad(A,fd);CHKERRQ(ierr);
  ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);

  ierr = MatCreateVecs(A,&x,&b);CHKERRQ(ierr);
  ierr = VecSetFromOptions(b);CHKERRQ(ierr);
  ierr = VecSet(b,1.0);CHKERRQ(ierr);

  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
  ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
  ierr = KSPSetUp(ksp);CHKERRQ(ierr);
  ierr = KSPSetUpOnBlocks(ksp);CHKERRQ(ierr);

  PetscPreLoadStage("KSPSolve");
  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);

  ierr = MatDestroy(&A);CHKERRQ(ierr);
  ierr = VecDestroy(&b);CHKERRQ(ierr);
  ierr = VecDestroy(&x);CHKERRQ(ierr);
  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
  PetscPreLoadEnd();
  ierr = PetscLogView_VecScatter(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

  ierr = PetscFinalize();
  return ierr;
}

int main(int argc, char **argv)
{
  PetscErrorCode ierr;
  Mat            A;
  KSP            ksp;
  PC             pc;
  IS             zero, one;
  MatNullSpace   nullsp;
  Vec            x, b;
  MPI_Comm       comm;

  PetscInitialize(&argc, &argv, NULL, NULL);

  comm = PETSC_COMM_WORLD;

  ierr = MatCreate(comm, &A);CHKERRQ(ierr);
  ierr = MatSetSizes(A, 4, 4, PETSC_DECIDE, PETSC_DECIDE);CHKERRQ(ierr);
  ierr = MatSetUp(A);CHKERRQ(ierr);
  ierr = MatSetFromOptions(A);CHKERRQ(ierr);
  ierr = MatCreateVecs(A, &x, &b);CHKERRQ(ierr);
  ierr = VecSet(x, 2.0);CHKERRQ(ierr);
  ierr = VecSet(b, 12.0);CHKERRQ(ierr);
  ierr = MatDiagonalSet(A, x, INSERT_VALUES);CHKERRQ(ierr);
  ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = ISCreateStride(comm, 2, 0, 1, &zero);CHKERRQ(ierr);
  ierr = ISCreateStride(comm, 2, 2, 1, &one);CHKERRQ(ierr);
  ierr = MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, &nullsp);CHKERRQ(ierr);
  ierr = PetscObjectCompose((PetscObject)zero, "nullspace",(PetscObject)nullsp);CHKERRQ(ierr);
  ierr = KSPCreate(comm, &ksp);CHKERRQ(ierr);
  ierr = KSPSetOperators(ksp, A, A);CHKERRQ(ierr);
  ierr = KSPSetUp(ksp);CHKERRQ(ierr);
  ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);
  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
  ierr = PCFieldSplitSetIS(pc, "0", zero);
  ierr = PCFieldSplitSetIS(pc, "1", one);
  ierr = KSPSolve(ksp, b, x);CHKERRQ(ierr);
  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
  ierr = MatNullSpaceDestroy(&nullsp);CHKERRQ(ierr);
  ierr = ISDestroy(&zero);CHKERRQ(ierr);
  ierr = ISDestroy(&one);CHKERRQ(ierr);
  ierr = MatDestroy(&A);CHKERRQ(ierr);
  ierr = VecDestroy(&x);CHKERRQ(ierr);
  ierr = VecDestroy(&b);CHKERRQ(ierr);

  PetscFinalize();
  return 0;
}

static PetscErrorCode KSPSetUp_SpecEst(KSP ksp)
{
  KSP_SpecEst    *spec = (KSP_SpecEst*)ksp->data;
  PetscErrorCode ierr;
  PetscBool      nonzero;

  PetscFunctionBegin;
  ierr = KSPSetPC(spec->kspest,ksp->pc);CHKERRQ(ierr);
  ierr = KSPSetPC(spec->kspcheap,ksp->pc);CHKERRQ(ierr);
  ierr = KSPGetInitialGuessNonzero(ksp,&nonzero);CHKERRQ(ierr);
  ierr = KSPSetInitialGuessNonzero(spec->kspest,nonzero);CHKERRQ(ierr);
  ierr = KSPSetInitialGuessNonzero(spec->kspcheap,nonzero);CHKERRQ(ierr);
  ierr = KSPSetComputeSingularValues(spec->kspest,PETSC_TRUE);CHKERRQ(ierr);
  ierr = KSPSetUp(spec->kspest);CHKERRQ(ierr);
  spec->current    = PETSC_FALSE;
  PetscFunctionReturn(0);
}