static PetscErrorCode KSPSolve_GCR(KSP ksp)
{
  KSP_GCR        *ctx = (KSP_GCR*)ksp->data;
  PetscErrorCode ierr;
  Mat            A, B;
  Vec            r,b,x;
  PetscReal      norm_r;

  PetscFunctionBegin;
  ierr = KSPGetOperators(ksp, &A, &B);CHKERRQ(ierr);
  x    = ksp->vec_sol;
  b    = ksp->vec_rhs;
  r    = ctx->R;

  /* compute initial residual */
  ierr = KSP_MatMult(ksp,A, x, r);CHKERRQ(ierr);
  ierr = VecAYPX(r, -1.0, b);CHKERRQ(ierr); /* r = b - A x  */
  ierr = VecNorm(r, NORM_2, &norm_r);CHKERRQ(ierr);
  KSPCheckNorm(ksp,norm_r);
  ksp->its    = 0;
  ksp->rnorm0 = norm_r;

  ierr = KSPLogResidualHistory(ksp,ksp->rnorm0);CHKERRQ(ierr);
  ierr = KSPMonitor(ksp,ksp->its,ksp->rnorm0);CHKERRQ(ierr);
  ierr = (*ksp->converged)(ksp,ksp->its,ksp->rnorm0,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
  if (ksp->reason) PetscFunctionReturn(0);

  do {
    ierr = KSPSolve_GCR_cycle(ksp);CHKERRQ(ierr);
    if (ksp->reason) break; /* catch case when convergence occurs inside the cycle */
  } while (ksp->its < ksp->max_it);CHKERRQ(ierr);

  if (ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
  PetscFunctionReturn(0);
}

static PetscErrorCode KSPSolve_CGLS(KSP ksp)
{
  PetscErrorCode ierr;
  KSP_CGLS       *cgls = (KSP_CGLS*)ksp->data;
  Mat            A;
  Vec            x,b,r,p,q,ss;
  PetscScalar    beta;
  PetscReal      alpha,gamma,oldgamma;
  PetscInt       maxiter_ls = 15;
  
  PetscFunctionBegin;
  ierr = KSPGetOperators(ksp,&A,NULL);CHKERRQ(ierr); /* Matrix of the system */
  
  /* vectors of length n, where system size is mxn */
  x  = ksp->vec_sol; /* Solution vector */
  p  = cgls->vwork_n[0];
  ss  = cgls->vwork_n[1];
  
  /* vectors of length m, where system size is mxn */
  b  = ksp->vec_rhs; /* Right-hand side vector */
  r  = cgls->vwork_m[0];
  q  = cgls->vwork_m[1];
  
  /* Minimization with the CGLS method */
  ksp->its = 0;
  ierr = MatMult(A,x,r);CHKERRQ(ierr);
  ierr = VecAYPX(r,-1,b);CHKERRQ(ierr);         /* r_0 = b - A * x_0  */
  ierr = MatMultTranspose(A,r,p);CHKERRQ(ierr); /* p_0 = A^T * r_0    */
  ierr = VecCopy(p,ss);CHKERRQ(ierr);           /* s_0 = p_0          */
  ierr = VecNorm(ss,NORM_2,&gamma);CHKERRQ(ierr);
  KSPCheckNorm(ksp,gamma);
  ksp->rnorm = gamma;
  ierr = (*ksp->converged)(ksp,ksp->its,ksp->rnorm,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
  gamma = gamma*gamma;                          /* gamma = norm2(s)^2 */

  do {
    ierr = MatMult(A,p,q);CHKERRQ(ierr);           /* q = A * p               */
    ierr = VecNorm(q,NORM_2,&alpha);CHKERRQ(ierr);
    KSPCheckNorm(ksp,alpha);
    alpha = alpha*alpha;                           /* alpha = norm2(q)^2      */
    alpha = gamma / alpha;                         /* alpha = gamma / alpha   */
    ierr = VecAXPY(x,alpha,p);CHKERRQ(ierr);       /* x += alpha * p          */
    ierr = VecAXPY(r,-alpha,q);CHKERRQ(ierr);      /* r -= alpha * q          */
    ierr = MatMultTranspose(A,r,ss);CHKERRQ(ierr); /* ss = A^T * r            */
    oldgamma = gamma;                              /* oldgamma = gamma        */
    ierr = VecNorm(ss,NORM_2,&gamma);CHKERRQ(ierr);
    KSPCheckNorm(ksp,gamma);
    ksp->its++;
    ksp->rnorm = gamma;
    ierr = KSPMonitor(ksp,ksp->its,ksp->rnorm);CHKERRQ(ierr);
    ierr = (*ksp->converged)(ksp,ksp->its,ksp->rnorm,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
    gamma = gamma*gamma;                           /* gamma = norm2(s)^2      */
    beta = gamma/oldgamma;                         /* beta = gamma / oldgamma */
    ierr = VecAYPX(p,beta,ss);CHKERRQ(ierr);       /* p = s + beta * p        */
  } while (ksp->its<ksp->max_it && !ksp->reason);
  
  if (ksp->its>=maxiter_ls && !ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
  PetscFunctionReturn(0);
}

static PetscErrorCode KSPAGMRESCycle(PetscInt *itcount,KSP ksp)
{
  KSP_AGMRES     *agmres = (KSP_AGMRES*)(ksp->data);
  PetscReal      res;
  PetscErrorCode ierr;
  PetscInt       KspSize = KSPSIZE;

  PetscFunctionBegin;
  /* check for the convergence */
  res = ksp->rnorm; /* Norm of the initial residual vector */
  if (!res) {
    if (itcount) *itcount = 0;
    ksp->reason = KSP_CONVERGED_ATOL;
    ierr        = PetscInfo(ksp,"Converged due to zero residual norm on entry\n");CHKERRQ(ierr);
    PetscFunctionReturn(0);
  }
  ierr = (*ksp->converged)(ksp,ksp->its,res,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
  /* Build the Krylov basis with Newton polynomials */
  ierr = KSPAGMRESBuildBasis(ksp);CHKERRQ(ierr);
  /* QR Factorize the basis with RODDEC */
  ierr = KSPAGMRESRoddec(ksp, KspSize+1);CHKERRQ(ierr);

  /* Recover a (partial) Hessenberg matrix for the Arnoldi-like relation */
  ierr = KSPAGMRESBuildHessenberg(ksp);CHKERRQ(ierr);
  /* Solve the least square problem and unwind the preconditioner */
  ierr = KSPAGMRESBuildSoln(ksp, KspSize);CHKERRQ(ierr);

  res        = ksp->rnorm;
  ksp->its  += KspSize;
  agmres->it = KspSize-1;
  /*  Test for the convergence */
  ierr = (*ksp->converged)(ksp,ksp->its,res,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
  ierr = KSPLogResidualHistory(ksp,res);CHKERRQ(ierr);
  ierr = KSPMonitor(ksp,ksp->its,res);CHKERRQ(ierr);


  *itcount = KspSize;
  PetscFunctionReturn(0);
}

PetscErrorCode KSPFGMRESCycle(PetscInt *itcount,KSP ksp)
{

  KSP_FGMRES     *fgmres = (KSP_FGMRES*)(ksp->data);
  PetscReal      res_norm;
  PetscReal      hapbnd,tt;
  PetscBool      hapend = PETSC_FALSE;  /* indicates happy breakdown ending */
  PetscErrorCode ierr;
  PetscInt       loc_it;                /* local count of # of dir. in Krylov space */
  PetscInt       max_k = fgmres->max_k; /* max # of directions Krylov space */
  Mat            Amat,Pmat;
  MatStructure   pflag;

  PetscFunctionBegin;
  /* Number of pseudo iterations since last restart is the number
     of prestart directions */
  loc_it = 0;

  /* note: (fgmres->it) is always set one less than (loc_it) It is used in
     KSPBUILDSolution_FGMRES, where it is passed to KSPFGMRESBuildSoln.
     Note that when KSPFGMRESBuildSoln is called from this function,
     (loc_it -1) is passed, so the two are equivalent */
  fgmres->it = (loc_it - 1);

  /* initial residual is in VEC_VV(0)  - compute its norm*/
  ierr = VecNorm(VEC_VV(0),NORM_2,&res_norm);CHKERRQ(ierr);

  /* first entry in right-hand-side of hessenberg system is just
     the initial residual norm */
  *RS(0) = res_norm;

  ksp->rnorm = res_norm;
  ierr       = KSPLogResidualHistory(ksp,res_norm);CHKERRQ(ierr);
  ierr       = KSPMonitor(ksp,ksp->its,res_norm);CHKERRQ(ierr);

  /* check for the convergence - maybe the current guess is good enough */
  ierr = (*ksp->converged)(ksp,ksp->its,res_norm,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
  if (ksp->reason) {
    if (itcount) *itcount = 0;
    PetscFunctionReturn(0);
  }

  /* scale VEC_VV (the initial residual) */
  ierr = VecScale(VEC_VV(0),1.0/res_norm);CHKERRQ(ierr);

  /* MAIN ITERATION LOOP BEGINNING*/
  /* keep iterating until we have converged OR generated the max number
     of directions OR reached the max number of iterations for the method */
  while (!ksp->reason && loc_it < max_k && ksp->its < ksp->max_it) {
    if (loc_it) {
      ierr = KSPLogResidualHistory(ksp,res_norm);CHKERRQ(ierr);
      ierr = KSPMonitor(ksp,ksp->its,res_norm);CHKERRQ(ierr);
    }
    fgmres->it = (loc_it - 1);

    /* see if more space is needed for work vectors */
    if (fgmres->vv_allocated <= loc_it + VEC_OFFSET + 1) {
      ierr = KSPFGMRESGetNewVectors(ksp,loc_it+1);CHKERRQ(ierr);
      /* (loc_it+1) is passed in as number of the first vector that should
         be allocated */
    }

    /* CHANGE THE PRECONDITIONER? */
    /* ModifyPC is the callback function that can be used to
       change the PC or its attributes before its applied */
    (*fgmres->modifypc)(ksp,ksp->its,loc_it,res_norm,fgmres->modifyctx);


    /* apply PRECONDITIONER to direction vector and store with
       preconditioned vectors in prevec */
    ierr = KSP_PCApply(ksp,VEC_VV(loc_it),PREVEC(loc_it));CHKERRQ(ierr);

    ierr = PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);CHKERRQ(ierr);
    /* Multiply preconditioned vector by operator - put in VEC_VV(loc_it+1) */
    ierr = MatMult(Amat,PREVEC(loc_it),VEC_VV(1+loc_it));CHKERRQ(ierr);


    /* update hessenberg matrix and do Gram-Schmidt - new direction is in
       VEC_VV(1+loc_it)*/
    ierr = (*fgmres->orthog)(ksp,loc_it);CHKERRQ(ierr);

    /* new entry in hessenburg is the 2-norm of our new direction */
    ierr = VecNorm(VEC_VV(loc_it+1),NORM_2,&tt);CHKERRQ(ierr);

    *HH(loc_it+1,loc_it)  = tt;
    *HES(loc_it+1,loc_it) = tt;

    /* Happy Breakdown Check */
    hapbnd = PetscAbsScalar((tt) / *RS(loc_it));
    /* RS(loc_it) contains the res_norm from the last iteration  */
    hapbnd = PetscMin(fgmres->haptol,hapbnd);
    if (tt > hapbnd) {
      /* scale new direction by its norm */
      ierr = VecScale(VEC_VV(loc_it+1),1.0/tt);CHKERRQ(ierr);
    } else {
      /* This happens when the solution is exactly reached. */
      /* So there is no new direction... */
      ierr   = VecSet(VEC_TEMP,0.0);CHKERRQ(ierr);     /* set VEC_TEMP to 0 */
      hapend = PETSC_TRUE;
    }
//.........这里部分代码省略.........

PetscErrorCode KSPSolve_QCG(KSP ksp)
{
/*
   Correpondence with documentation above:
      B = g = gradient,
      X = s = step
   Note:  This is not coded correctly for complex arithmetic!
 */

  KSP_QCG        *pcgP = (KSP_QCG*)ksp->data;
  Mat            Amat,Pmat;
  Vec            W,WA,WA2,R,P,ASP,BS,X,B;
  PetscScalar    scal,beta,rntrn,step;
  PetscReal      q1,q2,xnorm,step1,step2,rnrm,btx,xtax;
  PetscReal      ptasp,rtr,wtasp,bstp;
  PetscReal      dzero = 0.0,bsnrm;
  PetscErrorCode ierr;
  PetscInt       i,maxit;
  PC             pc = ksp->pc;
  PCSide         side;
  PetscBool      diagonalscale;

  PetscFunctionBegin;
  ierr = PCGetDiagonalScale(ksp->pc,&diagonalscale);CHKERRQ(ierr);
  if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
  if (ksp->transpose_solve) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Currently does not support transpose solve");

  ksp->its = 0;
  maxit    = ksp->max_it;
  WA       = ksp->work[0];
  R        = ksp->work[1];
  P        = ksp->work[2];
  ASP      = ksp->work[3];
  BS       = ksp->work[4];
  W        = ksp->work[5];
  WA2      = ksp->work[6];
  X        = ksp->vec_sol;
  B        = ksp->vec_rhs;

  if (pcgP->delta <= dzero) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Input error: delta <= 0");
  ierr = KSPGetPCSide(ksp,&side);CHKERRQ(ierr);
  if (side != PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Requires symmetric preconditioner!");

  /* Initialize variables */
  ierr = VecSet(W,0.0);CHKERRQ(ierr);  /* W = 0 */
  ierr = VecSet(X,0.0);CHKERRQ(ierr);  /* X = 0 */
  ierr = PCGetOperators(pc,&Amat,&Pmat);CHKERRQ(ierr);

  /* Compute:  BS = D^{-1} B */
  ierr = PCApplySymmetricLeft(pc,B,BS);CHKERRQ(ierr);

  ierr       = VecNorm(BS,NORM_2,&bsnrm);CHKERRQ(ierr);
  ierr       = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr);
  ksp->its   = 0;
  ksp->rnorm = bsnrm;
  ierr       = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr);
  ierr = KSPLogResidualHistory(ksp,bsnrm);CHKERRQ(ierr);
  ierr = KSPMonitor(ksp,0,bsnrm);CHKERRQ(ierr);
  ierr = (*ksp->converged)(ksp,0,bsnrm,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
  if (ksp->reason) PetscFunctionReturn(0);

  /* Compute the initial scaled direction and scaled residual */
  ierr = VecCopy(BS,R);CHKERRQ(ierr);
  ierr = VecScale(R,-1.0);CHKERRQ(ierr);
  ierr = VecCopy(R,P);CHKERRQ(ierr);
  ierr = VecDotRealPart(R,R,&rtr);CHKERRQ(ierr);

  for (i=0; i<=maxit; i++) {
    ierr = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr);
    ksp->its++;
    ierr = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr);

    /* Compute:  asp = D^{-T}*A*D^{-1}*p  */
    ierr = PCApplySymmetricRight(pc,P,WA);CHKERRQ(ierr);
    ierr = KSP_MatMult(ksp,Amat,WA,WA2);CHKERRQ(ierr);
    ierr = PCApplySymmetricLeft(pc,WA2,ASP);CHKERRQ(ierr);

    /* Check for negative curvature */
    ierr = VecDotRealPart(P,ASP,&ptasp);CHKERRQ(ierr);
    if (ptasp <= dzero) {

      /* Scaled negative curvature direction:  Compute a step so that
        ||w + step*p|| = delta and QS(w + step*p) is least */

      if (!i) {
        ierr = VecCopy(P,X);CHKERRQ(ierr);
        ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr);
        scal = pcgP->delta / xnorm;
        ierr = VecScale(X,scal);CHKERRQ(ierr);
      } else {
        /* Compute roots of quadratic */
        ierr = KSPQCGQuadraticRoots(W,P,pcgP->delta,&step1,&step2);CHKERRQ(ierr);
        ierr = VecDotRealPart(W,ASP,&wtasp);CHKERRQ(ierr);
        ierr = VecDotRealPart(BS,P,&bstp);CHKERRQ(ierr);
        ierr = VecCopy(W,X);CHKERRQ(ierr);
        q1   = step1*(bstp + wtasp + .5*step1*ptasp);
        q2   = step2*(bstp + wtasp + .5*step2*ptasp);
        if (q1 <= q2) {
          ierr = VecAXPY(X,step1,P);CHKERRQ(ierr);
        } else {
//.........这里部分代码省略.........

static PetscErrorCode KSPPGMRESCycle(PetscInt *itcount,KSP ksp)
{
  KSP_PGMRES     *pgmres = (KSP_PGMRES*)(ksp->data);
  PetscReal      res_norm,res,newnorm;
  PetscErrorCode ierr;
  PetscInt       it     = 0,j,k;
  PetscBool      hapend = PETSC_FALSE;

  PetscFunctionBegin;
  if (itcount) *itcount = 0;
  ierr   = VecNormalize(VEC_VV(0),&res_norm);CHKERRQ(ierr);
  res    = res_norm;
  *RS(0) = res_norm;

  /* check for the convergence */
  ierr       = PetscObjectAMSTakeAccess((PetscObject)ksp);CHKERRQ(ierr);
  ksp->rnorm = res;
  ierr       = PetscObjectAMSGrantAccess((PetscObject)ksp);CHKERRQ(ierr);
  pgmres->it = it-2;
  ierr = KSPLogResidualHistory(ksp,res);CHKERRQ(ierr);
  ierr = KSPMonitor(ksp,ksp->its,res);CHKERRQ(ierr);
  if (!res) {
    ksp->reason = KSP_CONVERGED_ATOL;
    ierr        = PetscInfo(ksp,"Converged due to zero residual norm on entry\n");CHKERRQ(ierr);
    PetscFunctionReturn(0);
  }

  ierr = (*ksp->converged)(ksp,ksp->its,res,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
  for (; !ksp->reason; it++) {
    Vec Zcur,Znext;
    if (pgmres->vv_allocated <= it + VEC_OFFSET + 1) {
      ierr = KSPGMRESGetNewVectors(ksp,it+1);CHKERRQ(ierr);
    }
    /* VEC_VV(it-1) is orthogonal, it will be normalized once the VecNorm arrives. */
    Zcur  = VEC_VV(it);         /* Zcur is not yet orthogonal, but the VecMDot to orthogonalize it has been started. */
    Znext = VEC_VV(it+1);       /* This iteration will compute Znext, update with a deferred correction once we know how
                                 * Zcur relates to the previous vectors, and start the reduction to orthogonalize it. */

    if (it < pgmres->max_k+1 && ksp->its+1 < PetscMax(2,ksp->max_it)) { /* We don't know whether what we have computed is enough, so apply the matrix. */
      ierr = KSP_PCApplyBAorAB(ksp,Zcur,Znext,VEC_TEMP_MATOP);CHKERRQ(ierr);
    }

    if (it > 1) {               /* Complete the pending reduction */
      ierr           = VecNormEnd(VEC_VV(it-1),NORM_2,&newnorm);CHKERRQ(ierr);
      *HH(it-1,it-2) = newnorm;
    }
    if (it > 0) {               /* Finish the reduction computing the latest column of H */
      ierr = VecMDotEnd(Zcur,it,&(VEC_VV(0)),HH(0,it-1));CHKERRQ(ierr);
    }

    if (it > 1) {
      /* normalize the base vector from two iterations ago, basis is complete up to here */
      ierr = VecScale(VEC_VV(it-1),1./ *HH(it-1,it-2));CHKERRQ(ierr);

      ierr       = KSPPGMRESUpdateHessenberg(ksp,it-2,&hapend,&res);CHKERRQ(ierr);
      pgmres->it = it-2;
      ksp->its++;
      ksp->rnorm = res;

      ierr = (*ksp->converged)(ksp,ksp->its,res,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
      if (it < pgmres->max_k+1 || ksp->reason || ksp->its == ksp->max_it) {  /* Monitor if we are done or still iterating, but not before a restart. */
        ierr = KSPLogResidualHistory(ksp,res);CHKERRQ(ierr);
        ierr = KSPMonitor(ksp,ksp->its,res);CHKERRQ(ierr);
      }
      if (ksp->reason) break;
      /* Catch error in happy breakdown and signal convergence and break from loop */
      if (hapend) {
        if (ksp->errorifnotconverged) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"You reached the happy break down, but convergence was not indicated. Residual norm = %G",res);
        else {
          ksp->reason = KSP_DIVERGED_BREAKDOWN;
          break;
        }
      }

      if (!(it < pgmres->max_k+1 && ksp->its < ksp->max_it)) break;

      /* The it-2 column of H was not scaled when we computed Zcur, apply correction */
      ierr = VecScale(Zcur,1./ *HH(it-1,it-2));CHKERRQ(ierr);
      /* And Znext computed in this iteration was computed using the under-scaled Zcur */
      ierr = VecScale(Znext,1./ *HH(it-1,it-2));CHKERRQ(ierr);

      /* In the previous iteration, we projected an unnormalized Zcur against the Krylov basis, so we need to fix the column of H resulting from that projection. */
      for (k=0; k<it; k++) *HH(k,it-1) /= *HH(it-1,it-2);
      /* When Zcur was projected against the Krylov basis, VV(it-1) was still not normalized, so fix that too. This
       * column is complete except for HH(it,it-1) which we won't know until the next iteration. */
      *HH(it-1,it-1) /= *HH(it-1,it-2);
    }

    if (it > 0) {
      PetscScalar *work;
      if (!pgmres->orthogwork) {ierr = PetscMalloc((pgmres->max_k + 2)*sizeof(PetscScalar),&pgmres->orthogwork);CHKERRQ(ierr);}
      work = pgmres->orthogwork;
      /* Apply correction computed by the VecMDot in the last iteration to Znext. The original form is
       *
       *   Znext -= sum_{j=0}^{i-1} Z[j+1] * H[j,i-1]
       *
       * where
       *
       *   Z[j] = sum_{k=0}^j V[k] * H[k,j-1]
       *
//.........这里部分代码省略.........

static PetscErrorCode  KSPSolve_CGS(KSP ksp)
{
  PetscErrorCode ierr;
  PetscInt       i;
  PetscScalar    rho,rhoold,a,s,b;
  Vec            X,B,V,P,R,RP,T,Q,U,AUQ;
  PetscReal      dp = 0.0;
  PetscBool      diagonalscale;

  PetscFunctionBegin;
  /* not sure what residual norm it does use, should use for right preconditioning */

  ierr = PCGetDiagonalScale(ksp->pc,&diagonalscale);CHKERRQ(ierr);
  if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);

  X   = ksp->vec_sol;
  B   = ksp->vec_rhs;
  R   = ksp->work[0];
  RP  = ksp->work[1];
  V   = ksp->work[2];
  T   = ksp->work[3];
  Q   = ksp->work[4];
  P   = ksp->work[5];
  U   = ksp->work[6];
  AUQ = V;

  /* Compute initial preconditioned residual */
  ierr = KSPInitialResidual(ksp,X,V,T,R,B);CHKERRQ(ierr);

  /* Test for nothing to do */
  ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr);
  if (ksp->normtype == KSP_NORM_NATURAL) dp *= dp;
  ierr       = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr);
  ksp->its   = 0;
  ksp->rnorm = dp;
  ierr       = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr);
  ierr = KSPLogResidualHistory(ksp,dp);CHKERRQ(ierr);
  ierr = KSPMonitor(ksp,0,dp);CHKERRQ(ierr);
  ierr = (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
  if (ksp->reason) PetscFunctionReturn(0);

  /* Make the initial Rp == R */
  ierr = VecCopy(R,RP);CHKERRQ(ierr);
  /*  added for Fidap */
  /* Penalize Startup - Isaac Hasbani Trick for CGS
     Since most initial conditions result in a mostly 0 residual,
     we change all the 0 values in the vector RP to the maximum.
  */
  if (ksp->normtype == KSP_NORM_NATURAL) {
    PetscReal   vr0max;
    PetscScalar *tmp_RP=0;
    PetscInt    numnp  =0, *max_pos=0;
    ierr = VecMax(RP, max_pos, &vr0max);CHKERRQ(ierr);
    ierr = VecGetArray(RP, &tmp_RP);CHKERRQ(ierr);
    ierr = VecGetLocalSize(RP, &numnp);CHKERRQ(ierr);
    for (i=0; i<numnp; i++) {
      if (tmp_RP[i] == 0.0) tmp_RP[i] = vr0max;
    }
    ierr = VecRestoreArray(RP, &tmp_RP);CHKERRQ(ierr);
  }
  /*  end of addition for Fidap */

  /* Set the initial conditions */
  ierr = VecDot(R,RP,&rhoold);CHKERRQ(ierr);        /* rhoold = (r,rp)      */
  ierr = VecCopy(R,U);CHKERRQ(ierr);
  ierr = VecCopy(R,P);CHKERRQ(ierr);
  ierr = KSP_PCApplyBAorAB(ksp,P,V,T);CHKERRQ(ierr);

  i = 0;
  do {

    ierr = VecDot(V,RP,&s);CHKERRQ(ierr);           /* s <- (v,rp)          */
    a    = rhoold / s;                               /* a <- rho / s         */
    ierr = VecWAXPY(Q,-a,V,U);CHKERRQ(ierr);      /* q <- u - a v         */
    ierr = VecWAXPY(T,1.0,U,Q);CHKERRQ(ierr);      /* t <- u + q           */
    ierr = VecAXPY(X,a,T);CHKERRQ(ierr);           /* x <- x + a (u + q)   */
    ierr = KSP_PCApplyBAorAB(ksp,T,AUQ,U);CHKERRQ(ierr);
    ierr = VecAXPY(R,-a,AUQ);CHKERRQ(ierr);       /* r <- r - a K (u + q) */
    ierr = VecDot(R,RP,&rho);CHKERRQ(ierr);         /* rho <- (r,rp)        */
    if (ksp->normtype == KSP_NORM_NATURAL) {
      dp = PetscAbsScalar(rho);
    } else {
      ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr);
    }

    ierr = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr);
    ksp->its++;
    ksp->rnorm = dp;
    ierr       = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr);
    ierr = KSPLogResidualHistory(ksp,dp);CHKERRQ(ierr);
    ierr = KSPMonitor(ksp,i+1,dp);CHKERRQ(ierr);
    ierr = (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
    if (ksp->reason) break;

    b      = rho / rhoold;                           /* b <- rho / rhoold    */
    ierr   = VecWAXPY(U,b,Q,R);CHKERRQ(ierr);       /* u <- r + b q         */
    ierr   = VecAXPY(Q,b,P);CHKERRQ(ierr);
    ierr   = VecWAXPY(P,b,Q,U);CHKERRQ(ierr);       /* p <- u + b(q + b p)  */
    ierr   = KSP_PCApplyBAorAB(ksp,P,V,Q);CHKERRQ(ierr);    /* v <- K p    */
    rhoold = rho;
//.........这里部分代码省略.........

PetscErrorCode KSPSolve_Chebyshev(KSP ksp)
{
  KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;
  PetscErrorCode ierr;
  PetscInt       k,kp1,km1,maxit,ktmp,i;
  PetscScalar    alpha,omegaprod,mu,omega,Gamma,c[3],scale;
  PetscReal      rnorm = 0.0;
  Vec            sol_orig,b,p[3],r;
  Mat            Amat,Pmat;
  PetscBool      diagonalscale;

  PetscFunctionBegin;
  ierr = PCGetDiagonalScale(ksp->pc,&diagonalscale);CHKERRQ(ierr);
  if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);

  if (cheb->kspest && !cheb->estimate_current) {
    PetscReal max=0.0,min=0.0;
    Vec       X,B;
    X = ksp->work[0];
    if (cheb->random) {
      B    = ksp->work[1];
      ierr = VecSetRandom(B,cheb->random);CHKERRQ(ierr);
    } else {
      B = ksp->vec_rhs;
    }
    ierr = KSPSolve(cheb->kspest,B,X);CHKERRQ(ierr);
    
    if (ksp->guess_zero) {
      ierr = VecZeroEntries(X);CHKERRQ(ierr);
    }
    ierr = KSPChebyshevComputeExtremeEigenvalues_Private(cheb->kspest,&min,&max);CHKERRQ(ierr);

    cheb->emin = cheb->tform[0]*min + cheb->tform[1]*max;
    cheb->emax = cheb->tform[2]*min + cheb->tform[3]*max;

    cheb->estimate_current = PETSC_TRUE;
  }

  ksp->its = 0;
  ierr     = PCGetOperators(ksp->pc,&Amat,&Pmat);CHKERRQ(ierr);
  maxit    = ksp->max_it;

  /* These three point to the three active solutions, we
     rotate these three at each solution update */
  km1      = 0; k = 1; kp1 = 2;
  sol_orig = ksp->vec_sol; /* ksp->vec_sol will be asigned to rotating vector p[k], thus save its address */
  b        = ksp->vec_rhs;
  p[km1]   = sol_orig;
  p[k]     = ksp->work[0];
  p[kp1]   = ksp->work[1];
  r        = ksp->work[2];

  /* use scale*B as our preconditioner */
  scale = 2.0/(cheb->emax + cheb->emin);

  /*   -alpha <=  scale*lambda(B^{-1}A) <= alpha   */
  alpha     = 1.0 - scale*(cheb->emin);
  Gamma     = 1.0;
  mu        = 1.0/alpha;
  omegaprod = 2.0/alpha;

  c[km1] = 1.0;
  c[k]   = mu;

  if (!ksp->guess_zero) {
    ierr = KSP_MatMult(ksp,Amat,p[km1],r);CHKERRQ(ierr);     /*  r = b - A*p[km1] */
    ierr = VecAYPX(r,-1.0,b);CHKERRQ(ierr);
  } else {
    ierr = VecCopy(b,r);CHKERRQ(ierr);
  }

  ierr = KSP_PCApply(ksp,r,p[k]);CHKERRQ(ierr);  /* p[k] = scale B^{-1}r + p[km1] */
  ierr = VecAYPX(p[k],scale,p[km1]);CHKERRQ(ierr);

  for (i=0; i<maxit; i++) {
    ierr = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr);

    ksp->its++;
    ierr   = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr);
    c[kp1] = 2.0*mu*c[k] - c[km1];
    omega  = omegaprod*c[k]/c[kp1];

    ierr = KSP_MatMult(ksp,Amat,p[k],r);CHKERRQ(ierr);          /*  r = b - Ap[k]    */
    ierr = VecAYPX(r,-1.0,b);CHKERRQ(ierr);
    ierr = KSP_PCApply(ksp,r,p[kp1]);CHKERRQ(ierr);             /*  p[kp1] = B^{-1}r  */
    ksp->vec_sol = p[k];

    /* calculate residual norm if requested */
    if (ksp->normtype != KSP_NORM_NONE || ksp->numbermonitors) {
      if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
        ierr = VecNorm(r,NORM_2,&rnorm);CHKERRQ(ierr);
      } else {
        ierr = VecNorm(p[kp1],NORM_2,&rnorm);CHKERRQ(ierr);
      }
      ierr         = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr);
      ksp->rnorm   = rnorm;
      ierr = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr);
      ierr = KSPLogResidualHistory(ksp,rnorm);CHKERRQ(ierr);
      ierr = KSPMonitor(ksp,i,rnorm);CHKERRQ(ierr);
      ierr = (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
//.........这里部分代码省略.........

static PetscErrorCode  KSPSolve_IBCGS(KSP ksp)
{
  PetscErrorCode ierr;
  PetscInt       i,N;
  PetscReal      rnorm,rnormin = 0.0;
#if defined(PETSC_HAVE_MPI_LONG_DOUBLE) && !defined(PETSC_USE_COMPLEX) && (defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL_DOUBLE))
  /* Because of possible instabilities in the algorithm (as indicated by different residual histories for the same problem 
     on the same number of processes  with different runs) we support computing the inner products using Intel's 80 bit arithematic
     rather than just 64 bit. Thus we copy our double precision values into long doubles (hoping this keeps the 16 extra bits)
     and tell MPI to do its ALlreduces with MPI_LONG_DOUBLE.

     Note for developers that does not effect the code. Intel's long double is implemented by storing the 80 bits of extended double
     precision into a 16 byte space (the rest of the space is ignored)  */
  long double    insums[7],outsums[7];
#else
  PetscScalar    insums[7],outsums[7];
#endif
  PetscScalar    sigman_2, sigman_1, sigman, pin_1, pin, phin_1, phin,tmp1,tmp2;
  PetscScalar    taun_1, taun, rhon, alphan_1, alphan, omegan_1, omegan;
  const PetscScalar *PETSC_RESTRICT r0, *PETSC_RESTRICT f0, *PETSC_RESTRICT qn, *PETSC_RESTRICT b, *PETSC_RESTRICT un;
  PetscScalar    *PETSC_RESTRICT rn, *PETSC_RESTRICT xn, *PETSC_RESTRICT vn, *PETSC_RESTRICT zn;
  /* the rest do not have to keep n_1 values */
  PetscScalar    kappan, thetan, etan, gamman, betan, deltan;
  const PetscScalar *PETSC_RESTRICT tn;
  PetscScalar    *PETSC_RESTRICT sn;
  Vec            R0,Rn,Xn,F0,Vn,Zn,Qn,Tn,Sn,B,Un;
  Mat            A;

  PetscFunctionBegin;
  if (!ksp->vec_rhs->petscnative) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"Only coded for PETSc vectors");

  ierr = PCGetOperators(ksp->pc,&A,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
  ierr = VecGetLocalSize(ksp->vec_sol,&N);CHKERRQ(ierr);
  Xn = ksp->vec_sol;ierr = VecGetArray(Xn_1,(PetscScalar**)&xn_1);CHKERRQ(ierr);ierr = VecRestoreArray(Xn_1,PETSC_NULL);CHKERRQ(ierr);
  B  = ksp->vec_rhs;ierr = VecGetArrayRead(B,(const PetscScalar**)&b);ierr = VecRestoreArrayRead(B,PETSC_NULL);CHKERRQ(ierr);
  R0 = ksp->work[0];ierr = VecGetArrayRead(R0,(const PetscScalar**)&r0);CHKERRQ(ierr);ierr = VecRestoreArrayRead(R0,PETSC_NULL);CHKERRQ(ierr);
  Rn = ksp->work[1];ierr = VecGetArray(Rn_1,(PetscScalar**)&rn_1);CHKERRQ(ierr);ierr = VecRestoreArray(Rn_1,PETSC_NULL);CHKERRQ(ierr);
  Un = ksp->work[2];ierr = VecGetArrayRead(Un_1,(const PetscScalar**)&un_1);CHKERRQ(ierr);ierr = VecRestoreArrayRead(Un_1,PETSC_NULL);CHKERRQ(ierr);
  F0 = ksp->work[3];ierr = VecGetArrayRead(F0,(const PetscScalar**)&f0);CHKERRQ(ierr);ierr = VecRestoreArrayRead(F0,PETSC_NULL);CHKERRQ(ierr);
  Vn = ksp->work[4];ierr = VecGetArray(Vn_1,(PetscScalar**)&vn_1);CHKERRQ(ierr);ierr = VecRestoreArray(Vn_1,PETSC_NULL);CHKERRQ(ierr);
  Zn = ksp->work[5];ierr = VecGetArray(Zn_1,(PetscScalar**)&zn_1);CHKERRQ(ierr);ierr = VecRestoreArray(Zn_1,PETSC_NULL);CHKERRQ(ierr);
  Qn = ksp->work[6];ierr = VecGetArrayRead(Qn_1,(const PetscScalar**)&qn_1);CHKERRQ(ierr);ierr = VecRestoreArrayRead(Qn_1,PETSC_NULL);CHKERRQ(ierr);
  Tn = ksp->work[7];ierr = VecGetArrayRead(Tn,(const PetscScalar**)&tn);CHKERRQ(ierr);ierr = VecRestoreArrayRead(Tn,PETSC_NULL);CHKERRQ(ierr);
  Sn = ksp->work[8];ierr = VecGetArrayRead(Sn,(const PetscScalar**)&sn);CHKERRQ(ierr);ierr = VecRestoreArrayRead(Sn,PETSC_NULL);CHKERRQ(ierr);

  /* r0 = rn_1 = b - A*xn_1; */
  /* ierr = KSP_PCApplyBAorAB(ksp,Xn_1,Rn_1,Tn);CHKERRQ(ierr);
     ierr = VecAYPX(Rn_1,-1.0,B);CHKERRQ(ierr); */
  ierr = KSPInitialResidual(ksp,Xn_1,Tn,Sn,Rn_1,B);CHKERRQ(ierr);

  ierr = VecNorm(Rn_1,NORM_2,&rnorm);CHKERRQ(ierr);
  ierr = KSPMonitor(ksp,0,rnorm);CHKERRQ(ierr);
  ierr = (*ksp->converged)(ksp,0,rnorm,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);   
  if (ksp->reason) PetscFunctionReturn(0);

  ierr = VecCopy(Rn_1,R0);CHKERRQ(ierr);

  /* un_1 = A*rn_1; */
  ierr = KSP_PCApplyBAorAB(ksp,Rn_1,Un_1,Tn);CHKERRQ(ierr);
  
  /* f0   = A'*rn_1; */
  if (ksp->pc_side == PC_RIGHT) { /* B' A' */
    ierr = MatMultTranspose(A,R0,Tn);CHKERRQ(ierr);
    ierr = PCApplyTranspose(ksp->pc,Tn,F0);CHKERRQ(ierr);
  } else if (ksp->pc_side == PC_LEFT) { /* A' B' */
    ierr = PCApplyTranspose(ksp->pc,R0,Tn);CHKERRQ(ierr);
    ierr = MatMultTranspose(A,Tn,F0);CHKERRQ(ierr);
  }

  /*qn_1 = vn_1 = zn_1 = 0.0; */
  ierr = VecSet(Qn_1,0.0);CHKERRQ(ierr);
  ierr = VecSet(Vn_1,0.0);CHKERRQ(ierr);
  ierr = VecSet(Zn_1,0.0);CHKERRQ(ierr);

  sigman_2 = pin_1 = taun_1 = 0.0;

  /* the paper says phin_1 should be initialized to zero, it is actually R0'R0 */
  ierr = VecDot(R0,R0,&phin_1);CHKERRQ(ierr); 

  /* sigman_1 = rn_1'un_1  */
  ierr = VecDot(R0,Un_1,&sigman_1);CHKERRQ(ierr); 

  alphan_1 = omegan_1 = 1.0;

  for (ksp->its = 1; ksp->its<ksp->max_it+1; ksp->its++) {
    rhon   = phin_1 - omegan_1*sigman_2 + omegan_1*alphan_1*pin_1;
    /*    if (rhon == 0.0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"rhon is zero, iteration %D",n); */
    if (ksp->its == 1) deltan = rhon;
    else deltan = rhon/taun_1;
    betan  = deltan/omegan_1;
    taun   = sigman_1 + betan*taun_1  - deltan*pin_1;
    if (taun == 0.0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"taun is zero, iteration %D",ksp->its);
    alphan = rhon/taun;
    ierr = PetscLogFlops(15.0);

    /*  
        zn = alphan*rn_1 + (alphan/alphan_1)betan*zn_1 - alphan*deltan*vn_1
        vn = un_1 + betan*vn_1 - deltan*qn_1
        sn = rn_1 - alphan*vn

//.........这里部分代码省略.........

PetscErrorCode  KSPSolve_BiCG(KSP ksp)
{
  PetscErrorCode ierr;
  PetscInt       i;
  PetscBool      diagonalscale;
  PetscScalar    dpi,a=1.0,beta,betaold=1.0,b,ma;
  PetscReal      dp;
  Vec            X,B,Zl,Zr,Rl,Rr,Pl,Pr;
  Mat            Amat,Pmat;

  PetscFunctionBegin;
  ierr = PCGetDiagonalScale(ksp->pc,&diagonalscale);CHKERRQ(ierr);
  if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);

  X  = ksp->vec_sol;
  B  = ksp->vec_rhs;
  Rl = ksp->work[0];
  Zl = ksp->work[1];
  Pl = ksp->work[2];
  Rr = ksp->work[3];
  Zr = ksp->work[4];
  Pr = ksp->work[5];

  ierr = PCGetOperators(ksp->pc,&Amat,&Pmat);CHKERRQ(ierr);

  if (!ksp->guess_zero) {
    ierr = KSP_MatMult(ksp,Amat,X,Rr);CHKERRQ(ierr);      /*   r <- b - Ax       */
    ierr = VecAYPX(Rr,-1.0,B);CHKERRQ(ierr);
  } else {
    ierr = VecCopy(B,Rr);CHKERRQ(ierr);           /*     r <- b (x is 0) */
  }
  ierr = VecCopy(Rr,Rl);CHKERRQ(ierr);
  ierr = KSP_PCApply(ksp,Rr,Zr);CHKERRQ(ierr);     /*     z <- Br         */
  ierr = VecConjugate(Rl);CHKERRQ(ierr);
  ierr = KSP_PCApplyTranspose(ksp,Rl,Zl);CHKERRQ(ierr);
  ierr = VecConjugate(Rl);CHKERRQ(ierr);
  ierr = VecConjugate(Zl);CHKERRQ(ierr);
  if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
    ierr = VecNorm(Zr,NORM_2,&dp);CHKERRQ(ierr);  /*    dp <- z'*z       */
  } else {
    ierr = VecNorm(Rr,NORM_2,&dp);CHKERRQ(ierr);  /*    dp <- r'*r       */
  }
  ierr       = KSPMonitor(ksp,0,dp);CHKERRQ(ierr);
  ierr       = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr);
  ksp->its   = 0;
  ksp->rnorm = dp;
  ierr       = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr);
  ierr = KSPLogResidualHistory(ksp,dp);CHKERRQ(ierr);
  ierr = (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
  if (ksp->reason) PetscFunctionReturn(0);

  i = 0;
  do {
    ierr = VecDot(Zr,Rl,&beta);CHKERRQ(ierr);       /*     beta <- r'z     */
    if (!i) {
      if (beta == 0.0) {
        ksp->reason = KSP_DIVERGED_BREAKDOWN_BICG;
        PetscFunctionReturn(0);
      }
      ierr = VecCopy(Zr,Pr);CHKERRQ(ierr);       /*     p <- z          */
      ierr = VecCopy(Zl,Pl);CHKERRQ(ierr);
    } else {
      b    = beta/betaold;
      ierr = VecAYPX(Pr,b,Zr);CHKERRQ(ierr);  /*     p <- z + b* p   */
      b    = PetscConj(b);
      ierr = VecAYPX(Pl,b,Zl);CHKERRQ(ierr);
    }
    betaold = beta;
    ierr    = KSP_MatMult(ksp,Amat,Pr,Zr);CHKERRQ(ierr); /*     z <- Kp         */
    ierr    = VecConjugate(Pl);CHKERRQ(ierr);
    ierr    = KSP_MatMultTranspose(ksp,Amat,Pl,Zl);CHKERRQ(ierr);
    ierr    = VecConjugate(Pl);CHKERRQ(ierr);
    ierr    = VecConjugate(Zl);CHKERRQ(ierr);
    ierr    = VecDot(Zr,Pl,&dpi);CHKERRQ(ierr);            /*     dpi <- z'p      */
    a       = beta/dpi;                           /*     a = beta/p'z    */
    ierr    = VecAXPY(X,a,Pr);CHKERRQ(ierr);    /*     x <- x + ap     */
    ma      = -a;
    ierr    = VecAXPY(Rr,ma,Zr);CHKERRQ(ierr);
    ma      = PetscConj(ma);
    ierr    = VecAXPY(Rl,ma,Zl);CHKERRQ(ierr);
    if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
      ierr = KSP_PCApply(ksp,Rr,Zr);CHKERRQ(ierr);  /*     z <- Br         */
      ierr = VecConjugate(Rl);CHKERRQ(ierr);
      ierr = KSP_PCApplyTranspose(ksp,Rl,Zl);CHKERRQ(ierr);
      ierr = VecConjugate(Rl);CHKERRQ(ierr);
      ierr = VecConjugate(Zl);CHKERRQ(ierr);
      ierr = VecNorm(Zr,NORM_2,&dp);CHKERRQ(ierr);  /*    dp <- z'*z       */
    } else {
      ierr = VecNorm(Rr,NORM_2,&dp);CHKERRQ(ierr);  /*    dp <- r'*r       */
    }
    ierr       = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr);
    ksp->its   = i+1;
    ksp->rnorm = dp;
    ierr       = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr);
    ierr = KSPLogResidualHistory(ksp,dp);CHKERRQ(ierr);
    ierr = KSPMonitor(ksp,i+1,dp);CHKERRQ(ierr);
    ierr = (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
    if (ksp->reason) break;
    if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
      ierr = KSP_PCApply(ksp,Rr,Zr);CHKERRQ(ierr);  /* z <- Br  */
//.........这里部分代码省略.........

static PetscErrorCode KSPPIPEFGMRESCycle(PetscInt *itcount,KSP ksp)
{
  KSP_PIPEFGMRES *pipefgmres = (KSP_PIPEFGMRES*)(ksp->data);
  PetscReal      res_norm;
  PetscReal      hapbnd,tt;
  PetscScalar    *hh,*hes,*lhh,shift = pipefgmres->shift;
  PetscBool      hapend = PETSC_FALSE;  /* indicates happy breakdown ending */
  PetscErrorCode ierr;
  PetscInt       loc_it;                /* local count of # of dir. in Krylov space */
  PetscInt       max_k = pipefgmres->max_k; /* max # of directions Krylov space */
  PetscInt       i,j,k;
  Mat            Amat,Pmat;
  Vec            Q,W; /* Pipelining vectors */
  Vec            *redux = pipefgmres->redux; /* workspace for single reduction */

  PetscFunctionBegin;
  if (itcount) *itcount = 0;

  /* Assign simpler names to these vectors, allocated as pipelining workspace */
  Q = VEC_Q;
  W = VEC_W;

  /* Allocate memory for orthogonalization work (freed in the GMRES Destroy routine)*/
  /* Note that we add an extra value here to allow for a single reduction */
  if (!pipefgmres->orthogwork) { ierr = PetscMalloc1(pipefgmres->max_k + 2 ,&pipefgmres->orthogwork);CHKERRQ(ierr);
  }
  lhh = pipefgmres->orthogwork;

  /* Number of pseudo iterations since last restart is the number
     of prestart directions */
  loc_it = 0;

  /* note: (pipefgmres->it) is always set one less than (loc_it) It is used in
     KSPBUILDSolution_PIPEFGMRES, where it is passed to KSPPIPEFGMRESBuildSoln.
     Note that when KSPPIPEFGMRESBuildSoln is called from this function,
     (loc_it -1) is passed, so the two are equivalent */
  pipefgmres->it = (loc_it - 1);

  /* initial residual is in VEC_VV(0)  - compute its norm*/
  ierr = VecNorm(VEC_VV(0),NORM_2,&res_norm);CHKERRQ(ierr);

  /* first entry in right-hand-side of hessenberg system is just
     the initial residual norm */
  *RS(0) = res_norm;

  ksp->rnorm = res_norm;
  ierr       = KSPLogResidualHistory(ksp,res_norm);CHKERRQ(ierr);
  ierr       = KSPMonitor(ksp,ksp->its,res_norm);CHKERRQ(ierr);

  /* check for the convergence - maybe the current guess is good enough */
  ierr = (*ksp->converged)(ksp,ksp->its,res_norm,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
  if (ksp->reason) {
    if (itcount) *itcount = 0;
    PetscFunctionReturn(0);
  }

  /* scale VEC_VV (the initial residual) */
  ierr = VecScale(VEC_VV(0),1.0/res_norm);CHKERRQ(ierr);

  /* Fill the pipeline */
  ierr = KSP_PCApply(ksp,VEC_VV(loc_it),PREVEC(loc_it));CHKERRQ(ierr);
  ierr = PCGetOperators(ksp->pc,&Amat,&Pmat);CHKERRQ(ierr);
  ierr = KSP_MatMult(ksp,Amat,PREVEC(loc_it),ZVEC(loc_it));CHKERRQ(ierr);
  ierr = VecAXPY(ZVEC(loc_it),-shift,VEC_VV(loc_it));CHKERRQ(ierr); /* Note shift */

  /* MAIN ITERATION LOOP BEGINNING*/
  /* keep iterating until we have converged OR generated the max number
     of directions OR reached the max number of iterations for the method */
  while (!ksp->reason && loc_it < max_k && ksp->its < ksp->max_it) {
    if (loc_it) {
      ierr = KSPLogResidualHistory(ksp,res_norm);CHKERRQ(ierr);
      ierr = KSPMonitor(ksp,ksp->its,res_norm);CHKERRQ(ierr);
    }
    pipefgmres->it = (loc_it - 1);

    /* see if more space is needed for work vectors */
    if (pipefgmres->vv_allocated <= loc_it + VEC_OFFSET + 1) {
      ierr = KSPPIPEFGMRESGetNewVectors(ksp,loc_it+1);CHKERRQ(ierr);
      /* (loc_it+1) is passed in as number of the first vector that should
         be allocated */
    }

    /* Note that these inner products are with "Z" now, so
       in particular, lhh[loc_it] is the 'barred' or 'shifted' value,
       not the value from the equivalent FGMRES run (even in exact arithmetic)
       That is, the H we need for the Arnoldi relation is different from the
       coefficients we use in the orthogonalization process,because of the shift */

    /* Do some local twiddling to allow for a single reduction */
    for(i=0;i<loc_it+1;i++){
      redux[i] = VEC_VV(i);
    }
    redux[loc_it+1] = ZVEC(loc_it);

    /* note the extra dot product which ends up in lh[loc_it+1], which computes ||z||^2 */
    ierr = VecMDotBegin(ZVEC(loc_it),loc_it+2,redux,lhh);CHKERRQ(ierr);

    /* Start the split reduction (This actually calls the MPI_Iallreduce, otherwise, the reduction is simply delayed until the "end" call)*/
    ierr = PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)ZVEC(loc_it)));CHKERRQ(ierr);

//.........这里部分代码省略.........

PetscErrorCode KSPSolve_TSIRM(KSP ksp)
{
    PetscErrorCode ierr;
    KSP_TSIRM      *tsirm = (KSP_TSIRM*)ksp->data;
    KSP            sub_ksp;
    PC             pc;
    Mat            AS;
    Vec            x,b;
    PetscScalar    *array;
    PetscReal      norm = 20;
    PetscInt       i,*ind_row,first_iteration = 1,its = 0,total = 0,col = 0;
    PetscInt       restart = 30;
    KSP            ksp_min;  /* KSP for minimization */
    PC             pc_min;    /* PC for minimization */

    PetscFunctionBegin;
    x = ksp->vec_sol; /* Solution vector        */
    b = ksp->vec_rhs; /* Right-hand side vector */

    /* Row indexes (these indexes are global) */
    ierr = PetscMalloc1(tsirm->Iend-tsirm->Istart,&ind_row);
    CHKERRQ(ierr);
    for (i=0; i<tsirm->Iend-tsirm->Istart; i++) ind_row[i] = i+tsirm->Istart;

    /* Inner solver */
    ierr = KSPGetPC(ksp,&pc);
    CHKERRQ(ierr);
    ierr = PCKSPGetKSP(pc,&sub_ksp);
    CHKERRQ(ierr);
    ierr = KSPSetTolerances(sub_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,restart);
    CHKERRQ(ierr);

    /* previously it seemed good but with SNES it seems not good... */
    ierr = KSP_MatMult(sub_ksp,tsirm->A,x,tsirm->r);
    CHKERRQ(ierr);
    ierr = VecAXPY(tsirm->r,-1,b);
    CHKERRQ(ierr);
    ierr = VecNorm(tsirm->r,NORM_2,&norm);
    CHKERRQ(ierr);
    ksp->its = 0;
    ierr = KSPConvergedDefault(ksp,ksp->its,norm,&ksp->reason,ksp->cnvP);
    CHKERRQ(ierr);
    ierr = KSPSetInitialGuessNonzero(sub_ksp,PETSC_TRUE);
    CHKERRQ(ierr);
    do {
        for (col=0; col<tsirm->size_ls && ksp->reason==0; col++) {
            /* Solve (inner iteration) */
            ierr = KSPSolve(sub_ksp,b,x);
            CHKERRQ(ierr);
            ierr = KSPGetIterationNumber(sub_ksp,&its);
            CHKERRQ(ierr);
            total += its;

            /* Build S^T */
            ierr = VecGetArray(x,&array);
            CHKERRQ(ierr);
            ierr = MatSetValues(tsirm->S,tsirm->Iend-tsirm->Istart,ind_row,1,&col,array,INSERT_VALUES);
            CHKERRQ(ierr);
            ierr = VecRestoreArray(x,&array);
            CHKERRQ(ierr);

            ierr = KSPGetResidualNorm(sub_ksp,&norm);
            CHKERRQ(ierr);
            ksp->rnorm = norm;
            ksp->its ++;
            ierr = KSPConvergedDefault(ksp,ksp->its,norm,&ksp->reason,ksp->cnvP);
            CHKERRQ(ierr);
            ierr = KSPMonitor(ksp,ksp->its,norm);
            CHKERRQ(ierr);
        }

        /* Minimization step */
        if (!ksp->reason) {
            ierr = MatAssemblyBegin(tsirm->S,MAT_FINAL_ASSEMBLY);
            CHKERRQ(ierr);
            ierr = MatAssemblyEnd(tsirm->S,MAT_FINAL_ASSEMBLY);
            CHKERRQ(ierr);
            if (first_iteration) {
                ierr = MatMatMult(tsirm->A,tsirm->S,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AS);
                CHKERRQ(ierr);
                first_iteration = 0;
            } else {
                ierr = MatMatMult(tsirm->A,tsirm->S,MAT_REUSE_MATRIX,PETSC_DEFAULT,&AS);
                CHKERRQ(ierr);
            }

            /* CGLS or LSQR method to minimize the residuals*/

            ierr = KSPCreate(PETSC_COMM_WORLD,&ksp_min);
            CHKERRQ(ierr);
            if (tsirm->cgls) {
                ierr = KSPSetType(ksp_min,KSPCGLS);
                CHKERRQ(ierr);
            } else {
                ierr = KSPSetType(ksp_min,KSPLSQR);
                CHKERRQ(ierr);
            }
            ierr = KSPSetOperators(ksp_min,AS,AS);
            CHKERRQ(ierr);
            ierr = KSPSetTolerances(ksp_min,tsirm->tol_ls,PETSC_DEFAULT,PETSC_DEFAULT,tsirm->maxiter_ls);
//.........这里部分代码省略.........

PetscErrorCode KSPLGMRESCycle(PetscInt *itcount,KSP ksp)
{
  KSP_LGMRES     *lgmres = (KSP_LGMRES*)(ksp->data);
  PetscReal      res_norm, res;
  PetscReal      hapbnd, tt;
  PetscScalar    tmp;
  PetscBool      hapend = PETSC_FALSE;  /* indicates happy breakdown ending */
  PetscErrorCode ierr;
  PetscInt       loc_it;                /* local count of # of dir. in Krylov space */
  PetscInt       max_k  = lgmres->max_k; /* max approx space size */
  PetscInt       max_it = ksp->max_it;  /* max # of overall iterations for the method */

  /* LGMRES_MOD - new variables*/
  PetscInt    aug_dim = lgmres->aug_dim;
  PetscInt    spot    = 0;
  PetscInt    order   = 0;
  PetscInt    it_arnoldi;                /* number of arnoldi steps to take */
  PetscInt    it_total;                  /* total number of its to take (=approx space size)*/
  PetscInt    ii, jj;
  PetscReal   tmp_norm;
  PetscScalar inv_tmp_norm;
  PetscScalar *avec;

  PetscFunctionBegin;
  /* Number of pseudo iterations since last restart is the number
     of prestart directions */