C++ PetscAbsScalar函数代码示例

static PetscErrorCode KSPSolve_PIPEFCG_cycle(KSP ksp)
{
  PetscErrorCode ierr;
  PetscInt       i,j,k,idx,kdx,mi;
  KSP_PIPEFCG    *pipefcg;
  PetscScalar    alpha=0.0,gamma,*betas,*dots;
  PetscReal      dp=0.0, delta,*eta,*etas;
  Vec            B,R,Z,X,Qcurr,W,ZETAcurr,M,N,Pcurr,Scurr,*redux;
  Mat            Amat,Pmat;

  PetscFunctionBegin;

  /* We have not checked these routines for use with complex numbers. The inner products
     are likely not defined correctly for that case */
#if (defined(PETSC_USE_COMPLEX) && !defined(PETSC_SKIP_COMPLEX))
  SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"PIPEFGMRES has not been implemented for use with complex scalars");
#endif

#define VecXDot(x,y,a)         (((pipefcg->type) == (KSP_CG_HERMITIAN)) ? VecDot       (x,y,a)   : VecTDot       (x,y,a))
#define VecXDotBegin(x,y,a)    (((pipefcg->type) == (KSP_CG_HERMITIAN)) ? VecDotBegin  (x,y,a)   : VecTDotBegin  (x,y,a))
#define VecXDotEnd(x,y,a)      (((pipefcg->type) == (KSP_CG_HERMITIAN)) ? VecDotEnd    (x,y,a)   : VecTDotEnd    (x,y,a))
#define VecMXDot(x,n,y,a)      (((pipefcg->type) == (KSP_CG_HERMITIAN)) ? VecMDot      (x,n,y,a) : VecMTDot      (x,n,y,a))
#define VecMXDotBegin(x,n,y,a) (((pipefcg->type) == (KSP_CG_HERMITIAN)) ? VecMDotBegin (x,n,y,a) : VecMTDotBegin (x,n,y,a))
#define VecMXDotEnd(x,n,y,a)   (((pipefcg->type) == (KSP_CG_HERMITIAN)) ? VecMDotEnd   (x,n,y,a) : VecMTDotEnd   (x,n,y,a))

  pipefcg       = (KSP_PIPEFCG*)ksp->data;
  X             = ksp->vec_sol;
  B             = ksp->vec_rhs;
  R             = ksp->work[0];
  Z             = ksp->work[1];
  W             = ksp->work[2];
  M             = ksp->work[3];
  N             = ksp->work[4];

  redux = pipefcg->redux;
  dots  = pipefcg->dots;
  etas  = pipefcg->etas;
  betas = dots;        /* dots takes the result of all dot products of which the betas are a subset */

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

  /* Compute cycle initial residual */
  ierr = KSP_MatMult(ksp,Amat,X,R);CHKERRQ(ierr);
  ierr = VecAYPX(R,-1.0,B);CHKERRQ(ierr);                   /* r <- b - Ax */
  ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr);                /* z <- Br     */

  Pcurr = pipefcg->Pvecs[0];
  Scurr = pipefcg->Svecs[0];
  Qcurr = pipefcg->Qvecs[0];
  ZETAcurr = pipefcg->ZETAvecs[0];
  ierr  = VecCopy(Z,Pcurr);CHKERRQ(ierr);
  ierr  = KSP_MatMult(ksp,Amat,Pcurr,Scurr);CHKERRQ(ierr);  /* S = Ap     */
  ierr  = VecCopy(Scurr,W);CHKERRQ(ierr);                   /* w = s = Az */

  /* Initial state of pipelining intermediates */
  redux[0] = R;
  redux[1] = W;
  ierr     = VecMXDotBegin(Z,2,redux,dots);CHKERRQ(ierr);
  ierr     = PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)Z));CHKERRQ(ierr); /* perform asynchronous reduction */
  ierr     = KSP_PCApply(ksp,W,M);CHKERRQ(ierr);            /* m = B(w) */
  ierr     = KSP_MatMult(ksp,Amat,M,N);CHKERRQ(ierr);       /* n = Am   */
  ierr     = VecCopy(M,Qcurr);CHKERRQ(ierr);                /* q = m    */
  ierr     = VecCopy(N,ZETAcurr);CHKERRQ(ierr);             /* zeta = n */
  ierr     = VecMXDotEnd(Z,2,redux,dots);CHKERRQ(ierr);
  gamma    = dots[0];
  delta    = PetscRealPart(dots[1]);
  etas[0]  = delta;
  alpha    = gamma/delta;

  i = 0;
  do {
    ksp->its++;

    /* Update X, R, Z, W */
    ierr = VecAXPY(X,+alpha,Pcurr);CHKERRQ(ierr);           /* x <- x + alpha * pi    */
    ierr = VecAXPY(R,-alpha,Scurr);CHKERRQ(ierr);           /* r <- r - alpha * si    */
    ierr = VecAXPY(Z,-alpha,Qcurr);CHKERRQ(ierr);           /* z <- z - alpha * qi    */
    ierr = VecAXPY(W,-alpha,ZETAcurr);CHKERRQ(ierr);        /* w <- w - alpha * zetai */

    /* Compute norm for convergence check */
    switch (ksp->normtype) {
      case KSP_NORM_PRECONDITIONED:
        ierr = VecNorm(Z,NORM_2,&dp);CHKERRQ(ierr);         /* dp <- sqrt(z'*z) = sqrt(e'*A'*B'*B*A*e) */
        break;
      case KSP_NORM_UNPRECONDITIONED:
        ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr);         /* dp <- sqrt(r'*r) = sqrt(e'*A'*A*e)      */
        break;
      case KSP_NORM_NATURAL:
        dp = PetscSqrtReal(PetscAbsScalar(gamma));          /* dp <- sqrt(r'*z) = sqrt(e'*A'*B*A*e)    */
        break;
      case KSP_NORM_NONE:
        dp = 0.0;
        break;
      default: SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"%s",KSPNormTypes[ksp->normtype]);
    }

    /* Check for convergence */
    ksp->rnorm = dp;
    KSPLogResidualHistory(ksp,dp);CHKERRQ(ierr);
    ierr = KSPMonitor(ksp,ksp->its,dp);CHKERRQ(ierr);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
          if (PetscRealPart(rval[k]) > 0) {
            a_pos += rval[k];
          } else {
            a_neg += rval[k];
          }
          ntotal++;
        } else diag = rval[k];
      }
      ierr = MatRestoreRow(lA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);

      /* ghosted all connections */
      if (gA) {
        ierr = MatGetRow(gA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
        for (k = 0; k < ncols; k++) {
          if (PetscRealPart(rval[k]) > 0.) {
            a_pos += PetscRealPart(rval[k]);
          } else {
            a_neg += PetscRealPart(rval[k]);
          }
          ntotal++;
        }
        ierr = MatRestoreRow(gA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
      }

      if (g_neg == 0.) {
        alpha = 0.;
      } else {
        alpha = -a_neg/g_neg;
      }

      if (g_pos == 0.) {
        diag += a_pos;
        beta = 0.;
      } else {
        beta = -a_pos/g_pos;
      }
      if (diag == 0.) {
        invdiag = 0.;
      } else invdiag = 1. / diag;
      /* on */
      ierr = MatGetRow(lG,i,&ncols,&icols,&vcols);CHKERRQ(ierr);
      idx = 0;
      for (j = 0;j < ncols;j++) {
        col = icols[j];
        if (lcid[col] >= 0 && vcols[j] != 0.) {
          row_f = i + fs;
          row_c = lcid[col];
          /* set the values for on-processor ones */
          if (PetscRealPart(vcols[j]) < 0.) {
            pij = vcols[j]*alpha*invdiag;
          } else {
            pij = vcols[j]*beta*invdiag;
          }
          if (PetscAbsScalar(pij) != 0.) {
            pvals[idx] = pij;
            pcols[idx] = row_c;
            idx++;
          }
        }
      }
      ierr = MatRestoreRow(lG,i,&ncols,&icols,&vcols);CHKERRQ(ierr);
      /* off */
      if (gG) {
        ierr = MatGetRow(gG,i,&ncols,&icols,&vcols);CHKERRQ(ierr);
        for (j = 0; j < ncols; j++) {
          col = icols[j];
          if (gcid[col] >= 0 && vcols[j] != 0.) {
            row_f = i + fs;
            row_c = gcid[col];
            /* set the values for on-processor ones */
            if (PetscRealPart(vcols[j]) < 0.) {
              pij = vcols[j]*alpha*invdiag;
            } else {
              pij = vcols[j]*beta*invdiag;
            }
            if (PetscAbsScalar(pij) != 0.) {
              pvals[idx] = pij;
              pcols[idx] = row_c;
              idx++;
            }
          }
        }
        ierr = MatRestoreRow(gG,i,&ncols,&icols,&vcols);CHKERRQ(ierr);
      }
      ierr = MatSetValues(*P,1,&row_f,idx,pcols,pvals,INSERT_VALUES);CHKERRQ(ierr);
    }
  }

  ierr = MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

  ierr = PetscFree(lsparse);CHKERRQ(ierr);
  ierr = PetscFree(gsparse);CHKERRQ(ierr);
  ierr = PetscFree(pcols);CHKERRQ(ierr);
  ierr = PetscFree(pvals);CHKERRQ(ierr);
  ierr = PetscFree(lcid);CHKERRQ(ierr);
  if (gG) {ierr = PetscFree(gcid);CHKERRQ(ierr);}

  PetscFunctionReturn(0);
}

PetscErrorCode  KSPSolve_CG(KSP ksp)
{
  PetscErrorCode ierr;
  PetscInt       i,stored_max_it,eigs;
  PetscScalar    dpi = 0.0,a = 1.0,beta,betaold = 1.0,b = 0,*e = 0,*d = 0,delta,dpiold;
  PetscReal      dp  = 0.0;
  Vec            X,B,Z,R,P,S,W;
  KSP_CG         *cg;
  Mat            Amat,Pmat;
  MatStructure   pflag;
  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);

  cg            = (KSP_CG*)ksp->data;
  eigs          = ksp->calc_sings;
  stored_max_it = ksp->max_it;
  X             = ksp->vec_sol;
  B             = ksp->vec_rhs;
  R             = ksp->work[0];
  Z             = ksp->work[1];
  P             = ksp->work[2];
  if (cg->singlereduction) {
    S = ksp->work[3];
    W = ksp->work[4];
  } else {
    S = 0;                      /* unused */
    W = Z;
  }

#define VecXDot(x,y,a) (((cg->type) == (KSP_CG_HERMITIAN)) ? VecDot(x,y,a) : VecTDot(x,y,a))

  if (eigs) {e = cg->e; d = cg->d; e[0] = 0.0; }
  ierr = PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);CHKERRQ(ierr);

  ksp->its = 0;
  if (!ksp->guess_zero) {
    ierr = KSP_MatMult(ksp,Amat,X,R);CHKERRQ(ierr);            /*     r <- b - Ax     */
    ierr = VecAYPX(R,-1.0,B);CHKERRQ(ierr);
  } else {
    ierr = VecCopy(B,R);CHKERRQ(ierr);                         /*     r <- b (x is 0) */
  }

  switch (ksp->normtype) {
  case KSP_NORM_PRECONDITIONED:
    ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr);                   /*     z <- Br         */
    ierr = VecNorm(Z,NORM_2,&dp);CHKERRQ(ierr);                /*    dp <- z'*z = e'*A'*B'*B*A'*e'     */
    break;
  case KSP_NORM_UNPRECONDITIONED:
    ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr);                /*    dp <- r'*r = e'*A'*A*e            */
    break;
  case KSP_NORM_NATURAL:
    ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr);                   /*     z <- Br         */
    if (cg->singlereduction) {
      ierr = KSP_MatMult(ksp,Amat,Z,S);CHKERRQ(ierr);
      ierr = VecXDot(Z,S,&delta);CHKERRQ(ierr);
    }
    ierr = VecXDot(Z,R,&beta);CHKERRQ(ierr);                     /*  beta <- z'*r       */
    if (PetscIsInfOrNanScalar(beta)) {
      if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf inner product");
      else {
        ksp->reason = KSP_DIVERGED_NANORINF;
        PetscFunctionReturn(0);
      }
    }
    dp = PetscSqrtReal(PetscAbsScalar(beta));                           /*    dp <- r'*z = r'*B*r = e'*A'*B*A*e */
    break;
  case KSP_NORM_NONE:
    dp = 0.0;
    break;
  default: SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"%s",KSPNormTypes[ksp->normtype]);
  }
  ierr       = KSPLogResidualHistory(ksp,dp);CHKERRQ(ierr);
  ierr       = KSPMonitor(ksp,0,dp);CHKERRQ(ierr);
  ksp->rnorm = dp;

  ierr = (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);      /* test for convergence */
  if (ksp->reason) PetscFunctionReturn(0);

  if (ksp->normtype != KSP_NORM_PRECONDITIONED && (ksp->normtype != KSP_NORM_NATURAL)) {
    ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr);                   /*     z <- Br         */
  }
  if (ksp->normtype != KSP_NORM_NATURAL) {
    if (cg->singlereduction) {
      ierr = KSP_MatMult(ksp,Amat,Z,S);CHKERRQ(ierr);
      ierr = VecXDot(Z,S,&delta);CHKERRQ(ierr);
    }
    ierr = VecXDot(Z,R,&beta);CHKERRQ(ierr);         /*  beta <- z'*r       */
    if (PetscIsInfOrNanScalar(beta)) {
      if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf inner product");
      else {
        ksp->reason = KSP_DIVERGED_NANORINF;
        PetscFunctionReturn(0);
      }
    }
  }

  i = 0;
//.........这里部分代码省略.........

static PetscErrorCode KSPLGMRESUpdateHessenberg(KSP ksp,PetscInt it,PetscBool  hapend,PetscReal *res)
{
    PetscScalar   *hh,*cc,*ss,tt;
    PetscInt      j;
    KSP_LGMRES    *lgmres = (KSP_LGMRES *)(ksp->data);

    PetscFunctionBegin;
    hh  = HH(0,it);  /* pointer to beginning of column to update - so
                      incrementing hh "steps down" the (it+1)th col of HH*/
    cc  = CC(0);     /* beginning of cosine rotations */
    ss  = SS(0);     /* beginning of sine rotations */

    /* Apply all the previously computed plane rotations to the new column
       of the Hessenberg matrix */
    /* Note: this uses the rotation [conj(c)  s ; -s   c], c= cos(theta), s= sin(theta) */

    for (j=1; j<=it; j++) {
        tt  = *hh;
        *hh = PetscConj(*cc) * tt + *ss * *(hh+1);
        hh++;
        *hh = *cc++ * *hh - (*ss++ * tt);
        /* hh, cc, and ss have all been incremented one by end of loop */
    }

    /*
      compute the new plane rotation, and apply it to:
       1) the right-hand-side of the Hessenberg system (GRS)
          note: it affects GRS(it) and GRS(it+1)
       2) the new column of the Hessenberg matrix
          note: it affects HH(it,it) which is currently pointed to
          by hh and HH(it+1, it) (*(hh+1))
      thus obtaining the updated value of the residual...
    */

    /* compute new plane rotation */

    if (!hapend) {
        tt        = PetscSqrtScalar(PetscConj(*hh) * *hh + PetscConj(*(hh+1)) * *(hh+1));
        if (tt == 0.0) {
            ksp->reason = KSP_DIVERGED_NULL;
            PetscFunctionReturn(0);
        }
        *cc       = *hh / tt;   /* new cosine value */
        *ss       = *(hh+1) / tt;  /* new sine value */

        /* apply to 1) and 2) */
        *GRS(it+1) = - (*ss * *GRS(it));
        *GRS(it)   = PetscConj(*cc) * *GRS(it);
        *hh        = PetscConj(*cc) * *hh + *ss * *(hh+1);

        /* residual is the last element (it+1) of right-hand side! */
        *res      = PetscAbsScalar(*GRS(it+1));

    } else {
        /* happy breakdown: HH(it+1, it) = 0, therfore we don't need to apply
                  another rotation matrix (so RH doesn't change).  The new residual is
                  always the new sine term times the residual from last time (GRS(it)),
                  but now the new sine rotation would be zero...so the residual should
                  be zero...so we will multiply "zero" by the last residual.  This might
                  not be exactly what we want to do here -could just return "zero". */

        *res = 0.0;
    }
    PetscFunctionReturn(0);
}

PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,Field **x,Field **f,void *ptr)
{
  AppCtx         *user = (AppCtx*)ptr;
  PetscErrorCode ierr;
  PetscInt       xints,xinte,yints,yinte,i,j;
  PetscReal      hx,hy,dhx,dhy,hxdhy,hydhx;
  PetscReal      grashof,prandtl,lid;
  PetscScalar    u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;

  PetscFunctionBeginUser;
  grashof = user->grashof;
  prandtl = user->prandtl;
  lid     = user->lidvelocity;

  /*
     Define mesh intervals ratios for uniform grid.

     Note: FD formulae below are normalized by multiplying through by
     local volume element (i.e. hx*hy) to obtain coefficients O(1) in two dimensions.


  */
  dhx   = (PetscReal)(info->mx-1);  dhy = (PetscReal)(info->my-1);
  hx    = 1.0/dhx;                   hy = 1.0/dhy;
  hxdhy = hx*dhy;                 hydhx = hy*dhx;

  xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;

  /* Test whether we are on the bottom edge of the global array */
  if (yints == 0) {
    j     = 0;
    yints = yints + 1;
    /* bottom edge */
    for (i=info->xs; i<info->xs+info->xm; i++) {
      f[j][i].u     = x[j][i].u;
      f[j][i].v     = x[j][i].v;
      f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
      f[j][i].temp  = x[j][i].temp-x[j+1][i].temp;
    }
  }

  /* Test whether we are on the top edge of the global array */
  if (yinte == info->my) {
    j     = info->my - 1;
    yinte = yinte - 1;
    /* top edge */
    for (i=info->xs; i<info->xs+info->xm; i++) {
      f[j][i].u     = x[j][i].u - lid;
      f[j][i].v     = x[j][i].v;
      f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
      f[j][i].temp  = x[j][i].temp-x[j-1][i].temp;
    }
  }

  /* Test whether we are on the left edge of the global array */
  if (xints == 0) {
    i     = 0;
    xints = xints + 1;
    /* left edge */
    for (j=info->ys; j<info->ys+info->ym; j++) {
      f[j][i].u     = x[j][i].u;
      f[j][i].v     = x[j][i].v;
      f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
      f[j][i].temp  = x[j][i].temp;
    }
  }

  /* Test whether we are on the right edge of the global array */
  if (xinte == info->mx) {
    i     = info->mx - 1;
    xinte = xinte - 1;
    /* right edge */
    for (j=info->ys; j<info->ys+info->ym; j++) {
      f[j][i].u     = x[j][i].u;
      f[j][i].v     = x[j][i].v;
      f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
      f[j][i].temp  = x[j][i].temp - (PetscReal)(grashof>0);
    }
  }

  /* Compute over the interior points */
  for (j=yints; j<yinte; j++) {
    for (i=xints; i<xinte; i++) {

      /*
       convective coefficients for upwinding
      */
      vx  = x[j][i].u; avx = PetscAbsScalar(vx);
      vxp = .5*(vx+avx); vxm = .5*(vx-avx);
      vy  = x[j][i].v; avy = PetscAbsScalar(vy);
      vyp = .5*(vy+avy); vym = .5*(vy-avy);

      /* U velocity */
      u         = x[j][i].u;
      uxx       = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
      uyy       = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
      f[j][i].u = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;

      /* V velocity */
      u         = x[j][i].v;
//.........这里部分代码省略.........

/*@C
   SNESComputeJacobianDefault - Computes the Jacobian using finite differences.

   Collective on SNES

   Input Parameters:
+  x1 - compute Jacobian at this point
-  ctx - application's function context, as set with SNESSetFunction()

   Output Parameters:
+  J - Jacobian matrix (not altered in this routine)
-  B - newly computed Jacobian matrix to use with preconditioner (generally the same as J)

   Options Database Key:
+  -snes_fd - Activates SNESComputeJacobianDefault()
.  -snes_test_err - Square root of function error tolerance, default square root of machine
                    epsilon (1.e-8 in double, 3.e-4 in single)
-  -mat_fd_type - Either wp or ds (see MATMFFD_WP or MATMFFD_DS)

   Notes:
   This routine is slow and expensive, and is not currently optimized
   to take advantage of sparsity in the problem.  Although
   SNESComputeJacobianDefault() is not recommended for general use
   in large-scale applications, It can be useful in checking the
   correctness of a user-provided Jacobian.

   An alternative routine that uses coloring to exploit matrix sparsity is
   SNESComputeJacobianDefaultColor().

   Level: intermediate

.keywords: SNES, finite differences, Jacobian

.seealso: SNESSetJacobian(), SNESComputeJacobianDefaultColor(), MatCreateSNESMF()
@*/
PetscErrorCode  SNESComputeJacobianDefault(SNES snes,Vec x1,Mat J,Mat B,void *ctx)
{
  Vec               j1a,j2a,x2;
  PetscErrorCode    ierr;
  PetscInt          i,N,start,end,j,value,root;
  PetscScalar       dx,*y,wscale;
  const PetscScalar *xx;
  PetscReal         amax,epsilon = PETSC_SQRT_MACHINE_EPSILON;
  PetscReal         dx_min = 1.e-16,dx_par = 1.e-1,unorm;
  MPI_Comm          comm;
  PetscErrorCode    (*eval_fct)(SNES,Vec,Vec)=0;
  PetscBool         assembled,use_wp = PETSC_TRUE,flg;
  const char        *list[2] = {"ds","wp"};
  PetscMPIInt       size;
  const PetscInt    *ranges;

  PetscFunctionBegin;
  ierr     = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_test_err",&epsilon,0);CHKERRQ(ierr);
  eval_fct = SNESComputeFunction;

  ierr = PetscObjectGetComm((PetscObject)x1,&comm);CHKERRQ(ierr);
  ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
  ierr = MatAssembled(B,&assembled);CHKERRQ(ierr);
  if (assembled) {
    ierr = MatZeroEntries(B);CHKERRQ(ierr);
  }
  if (!snes->nvwork) {
    snes->nvwork = 3;

    ierr = VecDuplicateVecs(x1,snes->nvwork,&snes->vwork);CHKERRQ(ierr);
    ierr = PetscLogObjectParents(snes,snes->nvwork,snes->vwork);CHKERRQ(ierr);
  }
  j1a = snes->vwork[0]; j2a = snes->vwork[1]; x2 = snes->vwork[2];

  ierr = VecGetSize(x1,&N);CHKERRQ(ierr);
  ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr);
  ierr = (*eval_fct)(snes,x1,j1a);CHKERRQ(ierr);

  ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"Differencing options","SNES");CHKERRQ(ierr);
  ierr = PetscOptionsEList("-mat_fd_type","Algorithm to compute difference parameter","SNESComputeJacobianDefault",list,2,"wp",&value,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsEnd();CHKERRQ(ierr);
  if (flg && !value) use_wp = PETSC_FALSE;

  if (use_wp) {
    ierr = VecNorm(x1,NORM_2,&unorm);CHKERRQ(ierr);
  }
  /* Compute Jacobian approximation, 1 column at a time.
      x1 = current iterate, j1a = F(x1)
      x2 = perturbed iterate, j2a = F(x2)
   */
  for (i=0; i<N; i++) {
    ierr = VecCopy(x1,x2);CHKERRQ(ierr);
    if (i>= start && i<end) {
      ierr = VecGetArrayRead(x1,&xx);CHKERRQ(ierr);
      if (use_wp) dx = PetscSqrtReal(1.0 + unorm);
      else        dx = xx[i-start];
      ierr = VecRestoreArrayRead(x1,&xx);CHKERRQ(ierr);
      if (PetscAbsScalar(dx) < dx_min) dx = (PetscRealPart(dx) < 0. ? -1. : 1.) * dx_par;
      dx    *= epsilon;
      wscale = 1.0/dx;
      ierr   = VecSetValues(x2,1,&i,&dx,ADD_VALUES);CHKERRQ(ierr);
    } else {
      wscale = 0.0;
    }
    ierr = VecAssemblyBegin(x2);CHKERRQ(ierr);
//.........这里部分代码省略.........

PetscErrorCode KSPSolve_FCG(KSP ksp)
{
  PetscErrorCode ierr;
  PetscInt       i,k,idx,mi;
  KSP_FCG        *fcg = (KSP_FCG*)ksp->data;
  PetscScalar    alpha=0.0,beta = 0.0,dpi,s;
  PetscReal      dp=0.0;
  Vec            B,R,Z,X,Pcurr,Ccurr;
  Mat            Amat,Pmat;
  PetscInt       eigs = ksp->calc_sings; /* Variables for eigen estimation - START*/
  PetscInt       stored_max_it = ksp->max_it;
  PetscScalar    alphaold = 0,betaold = 1.0,*e = 0,*d = 0;/* Variables for eigen estimation  - FINISH */

  PetscFunctionBegin;

#define VecXDot(x,y,a) (((fcg->type) == (KSP_CG_HERMITIAN)) ? VecDot(x,y,a) : VecTDot(x,y,a))
#define VecXMDot(a,b,c,d) (((fcg->type) == (KSP_CG_HERMITIAN)) ? VecMDot(a,b,c,d) : VecMTDot(a,b,c,d))

  X             = ksp->vec_sol;
  B             = ksp->vec_rhs;
  R             = ksp->work[0];
  Z             = ksp->work[1];

  ierr = PCGetOperators(ksp->pc,&Amat,&Pmat);CHKERRQ(ierr);
  if (eigs) {e = fcg->e; d = fcg->d; e[0] = 0.0; }
  /* Compute initial residual needed for convergence check*/
  ksp->its = 0;
  if (!ksp->guess_zero) {
    ierr = KSP_MatMult(ksp,Amat,X,R);CHKERRQ(ierr);
    ierr = VecAYPX(R,-1.0,B);CHKERRQ(ierr);                    /*   r <- b - Ax     */
  } else {
    ierr = VecCopy(B,R);CHKERRQ(ierr);                         /*   r <- b (x is 0) */
  }
  switch (ksp->normtype) {
    case KSP_NORM_PRECONDITIONED:
      ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr);               /*   z <- Br         */
      ierr = VecNorm(Z,NORM_2,&dp);CHKERRQ(ierr);              /*   dp <- dqrt(z'*z) = sqrt(e'*A'*B'*B*A*e)     */
      break;
    case KSP_NORM_UNPRECONDITIONED:
      ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr);              /*   dp <- sqrt(r'*r) = sqrt(e'*A'*A*e)     */
      break;
    case KSP_NORM_NATURAL:
      ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr);               /*   z <- Br         */
      ierr = VecXDot(R,Z,&s);CHKERRQ(ierr);
      dp = PetscSqrtReal(PetscAbsScalar(s));                   /*   dp <- sqrt(r'*z) = sqrt(e'*A'*B*A*e)  */
      break;
    case KSP_NORM_NONE:
      dp = 0.0;
      break;
    default: SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"%s",KSPNormTypes[ksp->normtype]);
  }

  /* Initial Convergence Check */
  ierr       = KSPLogResidualHistory(ksp,dp);CHKERRQ(ierr);
  ierr       = KSPMonitor(ksp,0,dp);CHKERRQ(ierr);
  ksp->rnorm = dp;
  if (ksp->normtype == KSP_NORM_NONE) {
    ierr = KSPConvergedSkip(ksp,0,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
  } else {
    ierr = (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
  }
  if (ksp->reason) PetscFunctionReturn(0);

  /* Apply PC if not already done for convergence check */
  if(ksp->normtype == KSP_NORM_UNPRECONDITIONED || ksp->normtype == KSP_NORM_NONE){
    ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr);               /*   z <- Br         */
  }

  i = 0;
  do {
    ksp->its = i+1;

    /*  If needbe, allocate a new chunk of vectors in P and C */
    ierr = KSPAllocateVectors_FCG(ksp,i+1,fcg->vecb);CHKERRQ(ierr);

    /* Note that we wrap around and start clobbering old vectors */
    idx = i % (fcg->mmax+1);
    Pcurr = fcg->Pvecs[idx];
    Ccurr = fcg->Cvecs[idx];

    /* Compute a new column of P (Currently does not support modified G-S or iterative refinement)*/
    switch(fcg->truncstrat){
      case KSP_FCG_TRUNC_TYPE_NOTAY :
        mi = PetscMax(1,i%(fcg->mmax+1));
        break;
      case KSP_FCG_TRUNC_TYPE_STANDARD :
        mi = fcg->mmax;
        break;
      default:
        SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unrecognized FCG Truncation Strategy");CHKERRQ(ierr);
    }
    ierr = VecCopy(Z,Pcurr);CHKERRQ(ierr);

    {
      PetscInt l,ndots;

      l = PetscMax(0,i-mi);
      ndots = i-l;
      if (ndots){
        PetscInt    j;
//.........这里部分代码省略.........

PETSC_EXTERN PetscErrorCode MatGetOrdering_AWBM(Mat A, MatOrderingType type, IS *permR, IS *permC)
{
  Vec *scalR, *scalC, scalRVec, scalCVec;
  scalR = &scalRVec; scalC = &scalCVec;

  /* EVERYTHING IS WRITTEN AS IF THE MATRIX WERE COLUMN-MAJOR */
  Mat_SeqAIJ      *aij = (Mat_SeqAIJ *) A->data;
  PetscInt         n   = A->rmap->n; /* Number of local columns */
  PetscInt         m   = A->cmap->n; /* Number of local rows */
  PetscInt        *match;            /* The row matched to each column, and inverse column permutation */
  PetscInt        *matchR;           /* The column matched to each row */
  PetscInt        *p;                /* The column permutation */
  const PetscInt  *ia  = aij->i;
  const PetscInt  *ja  = aij->j;
  const MatScalar *a   = aij->a;
  Vec              colMax;
  PetscScalar     *a_j, *sr, *sc;
  PetscReal       *weights /* c_ij */, *u /* u_i */, *v /* v_j */, eps = PETSC_SQRT_MACHINE_EPSILON;
  PetscInt         debug = 0, r, c, r1, c1;
  PetscErrorCode   ierr;

  PetscFunctionBegin;
  ierr = PetscOptionsGetInt(NULL, "-debug", &debug, NULL);CHKERRQ(ierr);
  ierr = MatGetVecs(A, NULL, &colMax);CHKERRQ(ierr);
  ierr = MatGetRowMaxAbs(A, colMax, NULL);CHKERRQ(ierr);
  ierr = PetscMalloc2(n, &match, m, &matchR);CHKERRQ(ierr);
  ierr = PetscMalloc1(n, &p);CHKERRQ(ierr);
  ierr = PetscCalloc3(m, &u, n, &v, ia[n], &weights);CHKERRQ(ierr);
  for (c = 0; c < n; ++c) match[c] = -1;
  /* Compute weights */
  ierr = VecGetArray(colMax, &a_j);CHKERRQ(ierr);
  for (c = 0; c < n; ++c) {
    for (r = ia[c]; r < ia[c+1]; ++r) {
      PetscReal ar = PetscAbsScalar(a[r]);

      if (ar == 0.0) weights[r] = PETSC_MAX_REAL;
      else           weights[r] = log(a_j[c]/ar);
    }
  }
  /* Compute local row weights */
  for (r = 0; r < m; ++r) u[r] = PETSC_MAX_REAL;
  for (c = 0; c < n; ++c) {
    for (r = ia[c]; r < ia[c+1]; ++r) {
      u[ja[r]] = PetscMin(u[ja[r]], weights[r]);
    }
  }
  /* Compute local column weights */
  for (c = 0; c < n; ++c) {
    v[c] = PETSC_MAX_REAL;
    for (r = ia[c]; r < ia[c+1]; ++r) {
      v[c] = PetscMin(v[c], weights[r] - u[ja[r]]);
    }
  }
  for (r = 0; r < m; ++r) matchR[r] = -1;
  /* Match columns */
  ierr = CheckUnmatched(n, match, matchR);CHKERRQ(ierr);
  for (c = 0; c < n; ++c) {
    /* if (match[c] >= 0) continue; */
    if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "Row %d\n  Weights:", c);CHKERRQ(ierr);}
    for (r = ia[c]; r < ia[c+1]; ++r) {
      PetscReal weight = weights[r] - u[ja[r]] - v[c];
      if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " %g", weight);CHKERRQ(ierr);}
      if ((weight <= eps) && (matchR[ja[r]] < 0)) {
        if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "Matched %d -- %d\n", c, ja[r]);CHKERRQ(ierr);}
        match[c]      = ja[r];
        matchR[ja[r]] = c;
        break;
      }
    }
    if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr);}
  }
  /* Deal with unmatched columns */
  ierr = CheckUnmatched(n, match, matchR);CHKERRQ(ierr);
  for (c = 0; c < n; ++c) {
    if (match[c] >= 0) continue;
    for (r = ia[c]; r < ia[c+1]; ++r) {
      PetscReal weight = weights[r] - u[ja[r]] - v[c];
      if (weight > eps) continue;
      /* \bar c_ij = 0 and (r, j1) \in M */
      c1 = matchR[ja[r]];
      for (r1 = ia[c1]; r1 < ia[c1+1]; ++r1) {
        PetscReal weight1 = weights[r1] - u[ja[r1]] - v[c1];
        if ((matchR[ja[r1]] < 0) && (weight1 <= eps)) {
          /* (r, c1) in M is replaced by (r, c) and (r1, c1) */
          if (debug) {
            ierr = PetscPrintf(PETSC_COMM_SELF, "Replaced match %d -- %d\n", c1, ja[r]);CHKERRQ(ierr);
            ierr = PetscPrintf(PETSC_COMM_SELF, "  Added  match %d -- %d\n", c,  ja[r]);CHKERRQ(ierr);
            ierr = PetscPrintf(PETSC_COMM_SELF, "  Added  match %d -- %d\n", c1, ja[r1]);CHKERRQ(ierr);
          }
          match[c]       = ja[r];
          matchR[ja[r]]  = c;
          match[c1]      = ja[r1];
          matchR[ja[r1]] = c1;
          break;
        }
      }
      if (match[c] >= 0) break;
    }
  }
  /* Allow matching with non-optimal rows */
//.........这里部分代码省略.........

/*
.  it - column of the Hessenberg that is complete, PGMRES is actually computing two columns ahead of this
 */
static PetscErrorCode KSPPGMRESUpdateHessenberg(KSP ksp,PetscInt it,PetscBool *hapend,PetscReal *res)
{
  PetscScalar    *hh,*cc,*ss,*rs;
  PetscInt       j;
  PetscReal      hapbnd;
  KSP_PGMRES     *pgmres = (KSP_PGMRES*)(ksp->data);
  PetscErrorCode ierr;

  PetscFunctionBegin;
  hh = HH(0,it);   /* pointer to beginning of column to update */
  cc = CC(0);      /* beginning of cosine rotations */
  ss = SS(0);      /* beginning of sine rotations */
  rs = RS(0);      /* right hand side of least squares system */

  /* The Hessenberg matrix is now correct through column it, save that form for possible spectral analysis */
  for (j=0; j<=it+1; j++) *HES(j,it) = hh[j];

  /* check for the happy breakdown */
  hapbnd = PetscMin(PetscAbsScalar(hh[it+1] / rs[it]),pgmres->haptol);
  if (PetscAbsScalar(hh[it+1]) < hapbnd) {
    ierr    = PetscInfo4(ksp,"Detected happy breakdown, current hapbnd = %14.12e H(%D,%D) = %14.12e\n",(double)hapbnd,it+1,it,(double)PetscAbsScalar(*HH(it+1,it)));CHKERRQ(ierr);
    *hapend = PETSC_TRUE;
  }

  /* Apply all the previously computed plane rotations to the new column
     of the Hessenberg matrix */
  /* Note: this uses the rotation [conj(c)  s ; -s   c], c= cos(theta), s= sin(theta),
     and some refs have [c   s ; -conj(s)  c] (don't be confused!) */

  for (j=0; j<it; j++) {
    PetscScalar hhj = hh[j];
    hh[j]   = PetscConj(cc[j])*hhj + ss[j]*hh[j+1];
    hh[j+1] =          -ss[j] *hhj + cc[j]*hh[j+1];
  }

  /*
    compute the new plane rotation, and apply it to:
     1) the right-hand-side of the Hessenberg system (RS)
        note: it affects RS(it) and RS(it+1)
     2) the new column of the Hessenberg matrix
        note: it affects HH(it,it) which is currently pointed to
        by hh and HH(it+1, it) (*(hh+1))
    thus obtaining the updated value of the residual...
  */

  /* compute new plane rotation */

  if (!*hapend) {
    PetscReal delta = PetscSqrtReal(PetscSqr(PetscAbsScalar(hh[it])) + PetscSqr(PetscAbsScalar(hh[it+1])));
    if (delta == 0.0) {
      ksp->reason = KSP_DIVERGED_NULL;
      PetscFunctionReturn(0);
    }

    cc[it] = hh[it] / delta;    /* new cosine value */
    ss[it] = hh[it+1] / delta;  /* new sine value */

    hh[it]   = PetscConj(cc[it])*hh[it] + ss[it]*hh[it+1];
    rs[it+1] = -ss[it]*rs[it];
    rs[it]   = PetscConj(cc[it])*rs[it];
    *res     = PetscAbsScalar(rs[it+1]);
  } else { /* happy breakdown: HH(it+1, it) = 0, therefore we don't need to apply
            another rotation matrix (so RH doesn't change).  The new residual is
            always the new sine term times the residual from last time (RS(it)),
            but now the new sine rotation would be zero...so the residual should
            be zero...so we will multiply "zero" by the last residual.  This might
            not be exactly what we want to do here -could just return "zero". */

    *res = 0.0;
  }
  PetscFunctionReturn(0);
}

static PetscErrorCode KSPSolve_TCQMR(KSP ksp)
{
  PetscReal      rnorm0,rnorm,dp1,Gamma;
  PetscScalar    theta,ep,cl1,sl1,cl,sl,sprod,tau_n1,f;
  PetscScalar    deltmp,rho,beta,eptmp,ta,s,c,tau_n,delta;
  PetscScalar    dp11,dp2,rhom1,alpha,tmp;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ksp->its = 0;

  ierr = KSPInitialResidual(ksp,x,u,v,r,b);CHKERRQ(ierr);
  ierr = VecNorm(r,NORM_2,&rnorm0);CHKERRQ(ierr);          /*  rnorm0 = ||r|| */

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

  ierr  = VecSet(um1,0.0);CHKERRQ(ierr);
  ierr  = VecCopy(r,u);CHKERRQ(ierr);
  rnorm = rnorm0;
  tmp   = 1.0/rnorm; ierr = VecScale(u,tmp);CHKERRQ(ierr);
  ierr  = VecSet(vm1,0.0);CHKERRQ(ierr);
  ierr  = VecCopy(u,v);CHKERRQ(ierr);
  ierr  = VecCopy(u,v0);CHKERRQ(ierr);
  ierr  = VecSet(pvec1,0.0);CHKERRQ(ierr);
  ierr  = VecSet(pvec2,0.0);CHKERRQ(ierr);
  ierr  = VecSet(p,0.0);CHKERRQ(ierr);
  theta = 0.0;
  ep    = 0.0;
  cl1   = 0.0;
  sl1   = 0.0;
  cl    = 0.0;
  sl    = 0.0;
  sprod = 1.0;
  tau_n1= rnorm0;
  f     = 1.0;
  Gamma = 1.0;
  rhom1 = 1.0;

  /*
   CALCULATE SQUARED LANCZOS  vectors
   */
  ierr = (*ksp->converged)(ksp,ksp->its,rnorm,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
  while (!ksp->reason) {
    ierr = KSPMonitor(ksp,ksp->its,rnorm);CHKERRQ(ierr);
    ksp->its++;

    ierr   = KSP_PCApplyBAorAB(ksp,u,y,vtmp);CHKERRQ(ierr); /* y = A*u */
    ierr   = VecDot(y,v0,&dp11);CHKERRQ(ierr);
    ierr   = VecDot(u,v0,&dp2);CHKERRQ(ierr);
    alpha  = dp11 / dp2;                          /* alpha = v0'*y/v0'*u */
    deltmp = alpha;
    ierr   = VecCopy(y,z);CHKERRQ(ierr);
    ierr   = VecAXPY(z,-alpha,u);CHKERRQ(ierr); /* z = y - alpha u */
    ierr   = VecDot(u,v0,&rho);CHKERRQ(ierr);
    beta   = rho / (f*rhom1);
    rhom1  = rho;
    ierr   = VecCopy(z,utmp);CHKERRQ(ierr);    /* up1 = (A-alpha*I)*
                                                (z-2*beta*p) + f*beta*
                                                beta*um1 */
    ierr   = VecAXPY(utmp,-2.0*beta,p);CHKERRQ(ierr);
    ierr   = KSP_PCApplyBAorAB(ksp,utmp,up1,vtmp);CHKERRQ(ierr);
    ierr   = VecAXPY(up1,-alpha,utmp);CHKERRQ(ierr);
    ierr   = VecAXPY(up1,f*beta*beta,um1);CHKERRQ(ierr);
    ierr   = VecNorm(up1,NORM_2,&dp1);CHKERRQ(ierr);
    f      = 1.0 / dp1;
    ierr   = VecScale(up1,f);CHKERRQ(ierr);
    ierr   = VecAYPX(p,-beta,z);CHKERRQ(ierr);   /* p = f*(z-beta*p) */
    ierr   = VecScale(p,f);CHKERRQ(ierr);
    ierr   = VecCopy(u,um1);CHKERRQ(ierr);
    ierr   = VecCopy(up1,u);CHKERRQ(ierr);
    beta   = beta/Gamma;
    eptmp  = beta;
    ierr   = KSP_PCApplyBAorAB(ksp,v,vp1,vtmp);CHKERRQ(ierr);
    ierr   = VecAXPY(vp1,-alpha,v);CHKERRQ(ierr);
    ierr   = VecAXPY(vp1,-beta,vm1);CHKERRQ(ierr);
    ierr   = VecNorm(vp1,NORM_2,&Gamma);CHKERRQ(ierr);
    ierr   = VecScale(vp1,1.0/Gamma);CHKERRQ(ierr);
    ierr   = VecCopy(v,vm1);CHKERRQ(ierr);
    ierr   = VecCopy(vp1,v);CHKERRQ(ierr);

    /*
       SOLVE  Ax = b
     */
    /* Apply last two Given's (Gl-1 and Gl) rotations to (beta,alpha,Gamma) */
    if (ksp->its > 2) {
      theta =  sl1*beta;
      eptmp = -cl1*beta;
    }
    if (ksp->its > 1) {
      ep     = -cl*eptmp + sl*alpha;
      deltmp = -sl*eptmp - cl*alpha;
    }
    if (PetscAbsReal(Gamma) > PetscAbsScalar(deltmp)) {
      ta = -deltmp / Gamma;
      s  = 1.0 / PetscSqrtScalar(1.0 + ta*ta);
      c  = s*ta;
    } else {
      ta = -Gamma/deltmp;
      c  = 1.0 / PetscSqrtScalar(1.0 + ta*ta);
//.........这里部分代码省略.........

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

        /*LGMRES_MOD: decide whether this is an arnoldi step or an aug step */
        if (loc_it < it_arnoldi) { /* Arnoldi */
            ierr = KSP_PCApplyBAorAB(ksp,VEC_VV(loc_it),VEC_VV(1+loc_it),VEC_TEMP_MATOP);
            CHKERRQ(ierr);
        } else { /*aug step */
            order = loc_it - it_arnoldi + 1; /* which aug step */
            for (ii=0; ii<aug_dim; ii++) {
                if (lgmres->aug_order[ii] == order) {
                    spot = ii;
                    break; /* must have this because there will be duplicates before aug_ct = aug_dim */
                }
            }

            ierr = VecCopy(A_AUGVEC(spot), VEC_VV(1+loc_it));
            CHKERRQ(ierr);
            /*note: an alternate implementation choice would be to only save the AUGVECS and
              not A_AUGVEC and then apply the PC here to the augvec */
        }

        /* update hessenberg matrix and do Gram-Schmidt - new direction is in
           VEC_VV(1+loc_it)*/
        ierr = (*lgmres->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;


        /* check for the happy breakdown */
        hapbnd  = PetscAbsScalar(tt / *GRS(loc_it));/* GRS(loc_it) contains the res_norm from the last iteration  */
        if (hapbnd > lgmres->haptol) hapbnd = lgmres->haptol;
        if (tt > hapbnd) {
            tmp = 1.0/tt;
            ierr = VecScale(VEC_VV(loc_it+1),tmp);
            CHKERRQ(ierr); /* scale new direction by its norm */
        } else {
            ierr = PetscInfo2(ksp,"Detected happy breakdown, current hapbnd = %G tt = %G\n",hapbnd,tt);
            CHKERRQ(ierr);
            hapend = PETSC_TRUE;
        }

        /* Now apply rotations to new col of hessenberg (and right side of system),
           calculate new rotation, and get new residual norm at the same time*/
        ierr = KSPLGMRESUpdateHessenberg(ksp,loc_it,hapend,&res);
        CHKERRQ(ierr);
        if (ksp->reason) break;

        loc_it++;
        lgmres->it  = (loc_it-1);  /* Add this here in case it has converged */

        ierr = PetscObjectTakeAccess(ksp);
        CHKERRQ(ierr);
        ksp->its++;
        ksp->rnorm = res;
        ierr = PetscObjectGrantAccess(ksp);
        CHKERRQ(ierr);

        ierr = (*ksp->converged)(ksp,ksp->its,res,&ksp->reason,ksp->cnvP);
        CHKERRQ(ierr);

        /* Catch error in happy breakdown and signal convergence and break from loop */
        if (hapend) {

PetscErrorCode  PCISApplyInvSchur(PC pc, Vec b, Vec x, Vec vec1_N, Vec vec2_N)
{
  PetscErrorCode ierr;
  PC_IS          *pcis = (PC_IS*)(pc->data);

  PetscFunctionBegin;
  /*
    Neumann solvers.
    Applying the inverse of the local Schur complement, i.e, solving a Neumann
    Problem with zero at the interior nodes of the RHS and extracting the interface
    part of the solution. inverse Schur complement is applied to b and the result
    is stored in x.
  */
  /* Setting the RHS vec1_N */
  ierr = VecSet(vec1_N,0.0);CHKERRQ(ierr);
  ierr = VecScatterBegin(pcis->N_to_B,b,vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
  ierr = VecScatterEnd  (pcis->N_to_B,b,vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
  /* Checking for consistency of the RHS */
  {
    PetscBool flg = PETSC_FALSE;
    ierr = PetscOptionsGetBool(NULL,"-pc_is_check_consistency",&flg,NULL);CHKERRQ(ierr);
    if (flg) {
      PetscScalar average;
      PetscViewer viewer;
      ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&viewer);CHKERRQ(ierr);

      ierr    = VecSum(vec1_N,&average);CHKERRQ(ierr);
      average = average / ((PetscReal)pcis->n);
      ierr    = PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);CHKERRQ(ierr);
      if (pcis->pure_neumann) {
        ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d is floating. Average = % 1.14e\n",PetscGlobalRank,PetscAbsScalar(average));CHKERRQ(ierr);
      } else {
        ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d is fixed.    Average = % 1.14e\n",PetscGlobalRank,PetscAbsScalar(average));CHKERRQ(ierr);
      }
      ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
      ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);CHKERRQ(ierr);
    }
  }
  /* Solving the system for vec2_N */
  ierr = KSPSolve(pcis->ksp_N,vec1_N,vec2_N);CHKERRQ(ierr);
  /* Extracting the local interface vector out of the solution */
  ierr = VecScatterBegin(pcis->N_to_B,vec2_N,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
  ierr = VecScatterEnd  (pcis->N_to_B,vec2_N,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

static PetscErrorCode KSPLGMRESBuildSoln(PetscScalar* nrs,Vec vguess,Vec vdest,KSP ksp,PetscInt it)
{
    PetscScalar    tt;
    PetscErrorCode ierr;
    PetscInt       ii,k,j;
    KSP_LGMRES     *lgmres = (KSP_LGMRES *)(ksp->data);
    /*LGMRES_MOD */
    PetscInt       it_arnoldi, it_aug;
    PetscInt       jj, spot = 0;

    PetscFunctionBegin;
    /* Solve for solution vector that minimizes the residual */

    /* If it is < 0, no lgmres steps have been performed */
    if (it < 0) {
        ierr = VecCopy(vguess,vdest);
        CHKERRQ(ierr); /* VecCopy() is smart, exists immediately if vguess == vdest */
        PetscFunctionReturn(0);
    }

    /* so (it+1) lgmres steps HAVE been performed */

    /* LGMRES_MOD - determine if we need to use augvecs for the soln  - do not assume that
       this is called after the total its allowed for an approx space */
    if (lgmres->approx_constant) {
        it_arnoldi = lgmres->max_k - lgmres->aug_ct;
    } else {
        it_arnoldi = lgmres->max_k - lgmres->aug_dim;
    }
    if (it_arnoldi >= it +1) {
        it_aug = 0;
        it_arnoldi = it+1;
    } else {
        it_aug = (it + 1) - it_arnoldi;
    }

    /* now it_arnoldi indicates the number of matvecs that took place */
    lgmres->matvecs += it_arnoldi;


    /* solve the upper triangular system - GRS is the right side and HH is
       the upper triangular matrix  - put soln in nrs */
    if (*HH(it,it) == 0.0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"HH(it,it) is identically zero; it = %D GRS(it) = %G",it,PetscAbsScalar(*GRS(it)));
    if (*HH(it,it) != 0.0) {
        nrs[it] = *GRS(it) / *HH(it,it);
    } else {
        nrs[it] = 0.0;
    }

    for (ii=1; ii<=it; ii++) {
        k   = it - ii;
        tt  = *GRS(k);
        for (j=k+1; j<=it; j++) tt  = tt - *HH(k,j) * nrs[j];
        nrs[k]   = tt / *HH(k,k);
    }

    /* Accumulate the correction to the soln of the preconditioned prob. in VEC_TEMP */
    ierr = VecSet(VEC_TEMP,0.0);
    CHKERRQ(ierr); /* set VEC_TEMP components to 0 */

    /*LGMRES_MOD - if augmenting has happened we need to form the solution
      using the augvecs */
    if (!it_aug) { /* all its are from arnoldi */
        ierr = VecMAXPY(VEC_TEMP,it+1,nrs,&VEC_VV(0));
        CHKERRQ(ierr);
    } else { /*use aug vecs */
        /*first do regular krylov directions */
        ierr = VecMAXPY(VEC_TEMP,it_arnoldi,nrs,&VEC_VV(0));
        CHKERRQ(ierr);
        /*now add augmented portions - add contribution of aug vectors one at a time*/


        for (ii=0; ii<it_aug; ii++) {
            for (jj=0; jj<lgmres->aug_dim; jj++) {
                if (lgmres->aug_order[jj] == (ii+1)) {
                    spot = jj;
                    break; /* must have this because there will be duplicates before aug_ct = aug_dim */
                }
            }
            ierr = VecAXPY(VEC_TEMP,nrs[it_arnoldi+ii],AUGVEC(spot));
            CHKERRQ(ierr);
        }
    }
    /* now VEC_TEMP is what we want to keep for augmenting purposes - grab before the
       preconditioner is "unwound" from right-precondtioning*/
    ierr = VecCopy(VEC_TEMP, AUG_TEMP);
    CHKERRQ(ierr);

    ierr = KSPUnwindPreconditioner(ksp,VEC_TEMP,VEC_TEMP_MATOP);
    CHKERRQ(ierr);

    /* add solution to previous solution */
    /* put updated solution into vdest.*/
    ierr = VecCopy(vguess,vdest);
    CHKERRQ(ierr);
    ierr = VecAXPY(vdest,1.0,VEC_TEMP);
    CHKERRQ(ierr);

    PetscFunctionReturn(0);
}

PetscErrorCode  KSPSolve_MINRES(KSP ksp)
{
  PetscErrorCode ierr;
  PetscInt       i;
  PetscScalar    alpha,beta,ibeta,betaold,eta,c=1.0,ceta,cold=1.0,coold,s=0.0,sold=0.0,soold;
  PetscScalar    rho0,rho1,irho1,rho2,mrho2,rho3,mrho3,dp = 0.0;
  PetscReal      np;
  Vec            X,B,R,Z,U,V,W,UOLD,VOLD,WOLD,WOOLD;
  Mat            Amat,Pmat;
  MatStructure   pflag;
  KSP_MINRES     *minres = (KSP_MINRES*)ksp->data;
  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);

  X     = ksp->vec_sol;
  B     = ksp->vec_rhs;
  R     = ksp->work[0];
  Z     = ksp->work[1];
  U     = ksp->work[2];
  V     = ksp->work[3];
  W     = ksp->work[4];
  UOLD  = ksp->work[5];
  VOLD  = ksp->work[6];
  WOLD  = ksp->work[7];
  WOOLD = ksp->work[8];

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

  ksp->its = 0;

  ierr = VecSet(UOLD,0.0);CHKERRQ(ierr);          /*     u_old  <-   0   */
  ierr = VecCopy(UOLD,VOLD);CHKERRQ(ierr);         /*     v_old  <-   0   */
  ierr = VecCopy(UOLD,W);CHKERRQ(ierr);            /*     w      <-   0   */
  ierr = VecCopy(UOLD,WOLD);CHKERRQ(ierr);         /*     w_old  <-   0   */

  if (!ksp->guess_zero) {
    ierr = KSP_MatMult(ksp,Amat,X,R);CHKERRQ(ierr); /*     r <- b - A*x    */
    ierr = VecAYPX(R,-1.0,B);CHKERRQ(ierr);
  } else {
    ierr = VecCopy(B,R);CHKERRQ(ierr);              /*     r <- b (x is 0) */
  }

  ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr); /*     z  <- B*r       */

  ierr = VecDot(R,Z,&dp);CHKERRQ(ierr);
  if (PetscRealPart(dp) < minres->haptol) {
    ierr = PetscInfo2(ksp,"Detected indefinite operator %G tolerance %G\n",PetscRealPart(dp),minres->haptol);CHKERRQ(ierr);
    ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
    PetscFunctionReturn(0);
  }

  dp   = PetscAbsScalar(dp);
  dp   = PetscSqrtScalar(dp);
  beta = dp;                                        /*  beta <- sqrt(r'*z  */
  eta  = beta;

  ierr  = VecCopy(R,V);CHKERRQ(ierr);
  ierr  = VecCopy(Z,U);CHKERRQ(ierr);
  ibeta = 1.0 / beta;
  ierr  = VecScale(V,ibeta);CHKERRQ(ierr);        /*    v <- r / beta     */
  ierr  = VecScale(U,ibeta);CHKERRQ(ierr);        /*    u <- z / beta     */

  ierr = VecNorm(Z,NORM_2,&np);CHKERRQ(ierr);      /*   np <- ||z||        */

  ierr       = KSPLogResidualHistory(ksp,np);CHKERRQ(ierr);
  ierr       = KSPMonitor(ksp,0,np);CHKERRQ(ierr);
  ksp->rnorm = np;
  ierr       = (*ksp->converged)(ksp,0,np,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); /* test for convergence */
  if (ksp->reason) PetscFunctionReturn(0);

  i = 0;
  do {
    ksp->its = i+1;

/*   Lanczos  */

    ierr = KSP_MatMult(ksp,Amat,U,R);CHKERRQ(ierr);   /*      r <- A*u   */
    ierr = VecDot(U,R,&alpha);CHKERRQ(ierr);          /*  alpha <- r'*u  */
    ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr); /*      z <- B*r   */

    ierr = VecAXPY(R,-alpha,V);CHKERRQ(ierr);     /*  r <- r - alpha v     */
    ierr = VecAXPY(Z,-alpha,U);CHKERRQ(ierr);     /*  z <- z - alpha u     */
    ierr = VecAXPY(R,-beta,VOLD);CHKERRQ(ierr);   /*  r <- r - beta v_old  */
    ierr = VecAXPY(Z,-beta,UOLD);CHKERRQ(ierr);   /*  z <- z - beta u_old  */

    betaold = beta;

    ierr = VecDot(R,Z,&dp);CHKERRQ(ierr);
    if ( PetscRealPart(dp) < minres->haptol) {
      ierr = PetscInfo2(ksp,"Detected indefinite operator %G tolerance %G\n",PetscRealPart(dp),minres->haptol);CHKERRQ(ierr);
      ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
      break;
    }

    dp   = PetscAbsScalar(dp);
    beta = PetscSqrtScalar(dp);                               /*  beta <- sqrt(r'*z)   */

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

int main(int argc,char **args)
{
  Mat             A,B;
  Vec             xx,s1,s2,yy;
  PetscErrorCode ierr;
  PetscInt        m=45,rows[2],cols[2],bs=1,i,row,col,*idx,M;
  PetscScalar     rval,vals1[4],vals2[4];
  PetscRandom     rdm;
  IS              is1,is2;
  PetscReal       s1norm,s2norm,rnorm,tol = 1.e-4;
  PetscBool       flg;
  MatFactorInfo   info;

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

  /* Test MatSetValues() and MatGetValues() */
  ierr = PetscOptionsGetInt(PETSC_NULL,"-mat_block_size",&bs,PETSC_NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetInt(PETSC_NULL,"-mat_size",&m,PETSC_NULL);CHKERRQ(ierr);
  M    = m*bs;
  ierr = MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,M,M,1,PETSC_NULL,&A);CHKERRQ(ierr);
  ierr = MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);
  ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,M,M,15,PETSC_NULL,&B);CHKERRQ(ierr);
  ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);
  ierr = PetscRandomCreate(PETSC_COMM_SELF,&rdm);CHKERRQ(ierr);
  ierr = PetscRandomSetFromOptions(rdm);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF,M,&xx);CHKERRQ(ierr);
  ierr = VecDuplicate(xx,&s1);CHKERRQ(ierr);
  ierr = VecDuplicate(xx,&s2);CHKERRQ(ierr);
  ierr = VecDuplicate(xx,&yy);CHKERRQ(ierr);

  /* For each row add atleast 15 elements */
  for (row=0; row<M; row++) {
    for (i=0; i<25*bs; i++) {
      ierr = PetscRandomGetValue(rdm,&rval);CHKERRQ(ierr);
      col  = (PetscInt)(PetscRealPart(rval)*M);
      ierr = MatSetValues(A,1,&row,1,&col,&rval,INSERT_VALUES);CHKERRQ(ierr);
      ierr = MatSetValues(B,1,&row,1,&col,&rval,INSERT_VALUES);CHKERRQ(ierr);
    }
  }

  /* Now set blocks of values */
  for (i=0; i<20*bs; i++) {
    ierr = PetscRandomGetValue(rdm,&rval);CHKERRQ(ierr);
    cols[0] = (PetscInt)(PetscRealPart(rval)*M);
    vals1[0] = rval;
    ierr = PetscRandomGetValue(rdm,&rval);CHKERRQ(ierr);
    cols[1] = (PetscInt)(PetscRealPart(rval)*M);
    vals1[1] = rval;
    ierr = PetscRandomGetValue(rdm,&rval);CHKERRQ(ierr);
    rows[0] = (PetscInt)(PetscRealPart(rval)*M);
    vals1[2] = rval;
    ierr = PetscRandomGetValue(rdm,&rval);CHKERRQ(ierr);
    rows[1] = (PetscInt)(PetscRealPart(rval)*M);
    vals1[3] = rval;
    ierr = MatSetValues(A,2,rows,2,cols,vals1,INSERT_VALUES);CHKERRQ(ierr);
    ierr = MatSetValues(B,2,rows,2,cols,vals1,INSERT_VALUES);CHKERRQ(ierr);
  }

  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

  /* Test MatNorm() */
  ierr = MatNorm(A,NORM_FROBENIUS,&s1norm);CHKERRQ(ierr);
  ierr = MatNorm(B,NORM_FROBENIUS,&s2norm);CHKERRQ(ierr);
  rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm;
  if ( rnorm>tol ) {
    ierr = PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_FROBENIUS()- NormA=%16.14e NormB=%16.14e bs = %D\n",s1norm,s2norm,bs);CHKERRQ(ierr);
  }
  ierr = MatNorm(A,NORM_INFINITY,&s1norm);CHKERRQ(ierr);
  ierr = MatNorm(B,NORM_INFINITY,&s2norm);CHKERRQ(ierr);
  rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm;
  if ( rnorm>tol ) {
    ierr = PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_INFINITY()- NormA=%16.14e NormB=%16.14e bs = %D\n",s1norm,s2norm,bs);CHKERRQ(ierr);
  }
  ierr = MatNorm(A,NORM_1,&s1norm);CHKERRQ(ierr);
  ierr = MatNorm(B,NORM_1,&s2norm);CHKERRQ(ierr);
  rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm;
  if ( rnorm>tol ) {
    ierr = PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_NORM_1()- NormA=%16.14e NormB=%16.14e bs = %D\n",s1norm,s2norm,bs);CHKERRQ(ierr);
  }

  /* MatShift() */
  rval = 10*s1norm;
  ierr = MatShift(A,rval);CHKERRQ(ierr);
  ierr = MatShift(B,rval);CHKERRQ(ierr);

  /* Test MatTranspose() */
  ierr = MatTranspose(A,MAT_REUSE_MATRIX,&A);CHKERRQ(ierr);
  ierr = MatTranspose(B,MAT_REUSE_MATRIX,&B);CHKERRQ(ierr);

  /* Now do MatGetValues()  */
  for (i=0; i<30; i++) {
    ierr = PetscRandomGetValue(rdm,&rval);CHKERRQ(ierr);
    cols[0] = (PetscInt)(PetscRealPart(rval)*M);
    ierr = PetscRandomGetValue(rdm,&rval);CHKERRQ(ierr);
    cols[1] = (PetscInt)(PetscRealPart(rval)*M);
    ierr = PetscRandomGetValue(rdm,&rval);CHKERRQ(ierr);
    rows[0] = (PetscInt)(PetscRealPart(rval)*M);
//.........这里部分代码省略.........