C++ VecCopy函数代码示例

PetscErrorCode SetInitialGuess(Vec X,Userctx *user)
{
  PetscErrorCode ierr;
  Vec            Xgen,Xnet;
  PetscScalar    *xgen,*xnet;
  PetscInt       i,idx=0;
  PetscScalar    Vr,Vi,IGr,IGi,Vm,Vm2;
  PetscScalar    Eqp,Edp,delta;
  PetscScalar    Efd,RF,VR; /* Exciter variables */
  PetscScalar    Id,Iq;  /* Generator dq axis currents */
  PetscScalar    theta,Vd,Vq,SE;

  PetscFunctionBegin;
  M[0] = 2*H[0]/w_s; M[1] = 2*H[1]/w_s; M[2] = 2*H[2]/w_s;
  D[0] = 0.1*M[0]; D[1] = 0.1*M[1]; D[2] = 0.1*M[2];

  ierr = DMCompositeGetLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr);

  /* Network subsystem initialization */
  ierr = VecCopy(user->V0,Xnet);CHKERRQ(ierr);

  /* Generator subsystem initialization */
  ierr = VecGetArray(Xgen,&xgen);CHKERRQ(ierr);
  ierr = VecGetArray(Xnet,&xnet);CHKERRQ(ierr);

  for (i=0; i < ngen; i++) {
    Vr  = xnet[2*gbus[i]]; /* Real part of generator terminal voltage */
    Vi  = xnet[2*gbus[i]+1]; /* Imaginary part of the generator terminal voltage */
    Vm  = PetscSqrtScalar(Vr*Vr + Vi*Vi); Vm2 = Vm*Vm;
    IGr = (Vr*PG[i] + Vi*QG[i])/Vm2;
    IGi = (Vi*PG[i] - Vr*QG[i])/Vm2;

    delta = atan2(Vi+Xq[i]*IGr,Vr-Xq[i]*IGi); /* Machine angle */

    theta = PETSC_PI/2.0 - delta;

    Id = IGr*PetscCosScalar(theta) - IGi*PetscSinScalar(theta); /* d-axis stator current */
    Iq = IGr*PetscSinScalar(theta) + IGi*PetscCosScalar(theta); /* q-axis stator current */

    Vd = Vr*PetscCosScalar(theta) - Vi*PetscSinScalar(theta);
    Vq = Vr*PetscSinScalar(theta) + Vi*PetscCosScalar(theta);

    Edp = Vd + Rs[i]*Id - Xqp[i]*Iq; /* d-axis transient EMF */
    Eqp = Vq + Rs[i]*Iq + Xdp[i]*Id; /* q-axis transient EMF */

    TM[i] = PG[i];

    /* The generator variables are ordered as [Eqp,Edp,delta,w,Id,Iq] */
    xgen[idx]   = Eqp;
    xgen[idx+1] = Edp;
    xgen[idx+2] = delta;
    xgen[idx+3] = w_s;

    idx = idx + 4;

    xgen[idx]   = Id;
    xgen[idx+1] = Iq;

    idx = idx + 2;

    /* Exciter */
    Efd = Eqp + (Xd[i] - Xdp[i])*Id;
    SE  = k1[i]*PetscExpScalar(k2[i]*Efd);
    VR  =  KE[i]*Efd + SE;
    RF  =  KF[i]*Efd/TF[i];

    xgen[idx]   = Efd;
    xgen[idx+1] = RF;
    xgen[idx+2] = VR;

    Vref[i] = Vm + (VR/KA[i]);

    idx = idx + 3;
  }

  ierr = VecRestoreArray(Xgen,&xgen);CHKERRQ(ierr);
  ierr = VecRestoreArray(Xnet,&xnet);CHKERRQ(ierr);

  /* ierr = VecView(Xgen,0);CHKERRQ(ierr); */
  ierr = DMCompositeGather(user->dmpgrid,X,INSERT_VALUES,Xgen,Xnet);CHKERRQ(ierr);
  ierr = DMCompositeRestoreLocalVectors(user->dmpgrid,&Xgen,&Xnet);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode SNESQNApply_Broyden(SNES snes,PetscInt it,Vec Y,Vec X,Vec Xold,Vec D)
{
  PetscErrorCode     ierr;
  SNES_QN            *qn = (SNES_QN*)snes->data;
  Vec                W   = snes->work[3];
  Vec                *U  = qn->U;
  PetscInt           m = qn->m;
  PetscInt           k,i,j,lits,l = m;
  PetscReal          unorm,a,b;
  PetscReal          *lambda=qn->lambda;
  PetscScalar        gdot;
  PetscReal          udot;

  PetscFunctionBegin;
  if (it < m) l = it;
  if (it > 0) {
    k = (it-1)%l;
    ierr = SNESLineSearchGetLambda(snes->linesearch,&lambda[k]);CHKERRQ(ierr);
    ierr = VecCopy(Xold,U[k]);CHKERRQ(ierr);
    ierr = VecAXPY(U[k],-1.0,X);CHKERRQ(ierr);
    if (qn->monitor) {
      ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
      ierr = PetscViewerASCIIPrintf(qn->monitor, "scaling vector %d of %d by lambda: %14.12e \n",k,l,lambda[k]);CHKERRQ(ierr);
      ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
    }
  }
  if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
    ierr = KSPSolve(snes->ksp,D,W);CHKERRQ(ierr);
    SNESCheckKSPSolve(snes);
    ierr              = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr);
    snes->linear_its += lits;
    ierr              = VecCopy(W,Y);CHKERRQ(ierr);
  } else {
    ierr = VecCopy(D,Y);CHKERRQ(ierr);
    ierr = VecScale(Y,qn->scaling);CHKERRQ(ierr);
  }

  /* inward recursion starting at the first update and working forward */
  for (i = 0; i < l-1; i++) {
    j = (it+i-l)%l;
    k = (it+i-l+1)%l;
    ierr = VecNorm(U[j],NORM_2,&unorm);CHKERRQ(ierr);
    ierr = VecDot(U[j],Y,&gdot);CHKERRQ(ierr);
    unorm *= unorm;
    udot = PetscRealPart(gdot);
    a = (lambda[j]/lambda[k]);
    b = -(1.-lambda[j]);
    a *= udot/unorm;
    b *= udot/unorm;
    ierr = VecAXPBYPCZ(Y,a,b,1.,U[k],U[j]);CHKERRQ(ierr);

    if (qn->monitor) {
      ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
      ierr = PetscViewerASCIIPrintf(qn->monitor, "using vector %d and %d, gdot: %14.12e\n",k,j,PetscRealPart(gdot));CHKERRQ(ierr);
      ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
    }
  }
  if (l > 0) {
    k = (it-1)%l;
    ierr = VecDot(U[k],Y,&gdot);CHKERRQ(ierr);
    ierr = VecNorm(U[k],NORM_2,&unorm);CHKERRQ(ierr);
    unorm *= unorm;
    udot = PetscRealPart(gdot);
    a = unorm/(unorm-lambda[k]*udot);
    b = -(1.-lambda[k])*udot/(unorm-lambda[k]*udot);
    if (qn->monitor) {
      ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
      ierr = PetscViewerASCIIPrintf(qn->monitor, "using vector %d: a: %14.12e b: %14.12e \n",k,a,b);CHKERRQ(ierr);
      ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
    }
    ierr = VecAXPBY(Y,b,a,U[k]);CHKERRQ(ierr);
  }
  l = m;
  if (it+1<m)l=it+1;
  k = it%l;
  if (qn->monitor) {
    ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
    ierr = PetscViewerASCIIPrintf(qn->monitor, "setting vector %d of %d\n",k,l);CHKERRQ(ierr);
    ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

PetscErrorCode SNESSolve_NEWTONLS(SNES snes)
{
  PetscErrorCode      ierr;
  PetscInt            maxits,i,lits;
  PetscBool           lssucceed;
  MatStructure        flg = DIFFERENT_NONZERO_PATTERN;
  PetscReal           fnorm,gnorm,xnorm,ynorm;
  Vec                 Y,X,F,G,W,FPC;
  KSPConvergedReason  kspreason;
  PetscBool           domainerror;
  SNESLineSearch      linesearch;
  SNESConvergedReason reason;

  PetscFunctionBegin;
  snes->numFailures            = 0;
  snes->numLinearSolveFailures = 0;
  snes->reason                 = SNES_CONVERGED_ITERATING;

  maxits        = snes->max_its;        /* maximum number of iterations */
  X             = snes->vec_sol;        /* solution vector */
  F             = snes->vec_func;       /* residual vector */
  Y             = snes->vec_sol_update; /* newton step */
  G             = snes->work[0];
  W             = snes->work[1];

  ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr);
  snes->iter = 0;
  snes->norm = 0.0;
  ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr);
  ierr = SNESGetSNESLineSearch(snes, &linesearch);CHKERRQ(ierr);
  if (!snes->vec_func_init_set) {
    ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr);
    ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr);
    if (domainerror) {
      snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
      PetscFunctionReturn(0);
    }
  } else {
    snes->vec_func_init_set = PETSC_FALSE;
  }
  if (!snes->norm_init_set) {
    ierr = VecNormBegin(F,NORM_2,&fnorm);CHKERRQ(ierr);        /* fnorm <- ||F||  */
    ierr = VecNormEnd(F,NORM_2,&fnorm);CHKERRQ(ierr);
    if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
  } else {
    fnorm = snes->norm_init;
    snes->norm_init_set = PETSC_FALSE;
  }
  ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr);
  snes->norm = fnorm;
  ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr);
  SNESLogConvHistory(snes,fnorm,0);
  ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr);

  /* set parameter for default relative tolerance convergence test */
  snes->ttol = fnorm*snes->rtol;
  /* test convergence */
  ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
  if (snes->reason) PetscFunctionReturn(0);

  for (i=0; i<maxits; i++) {

    /* Call general purpose update function */
    if (snes->ops->update) {
      ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr);
    }

    /* apply the nonlinear preconditioner if it's right preconditioned */
    if (snes->pc && snes->pcside == PC_RIGHT) {
      ierr = SNESSetInitialFunction(snes->pc, F);CHKERRQ(ierr);
      ierr = SNESSetInitialFunctionNorm(snes->pc, fnorm);CHKERRQ(ierr);
      ierr = SNESSolve(snes->pc, snes->vec_rhs, X);CHKERRQ(ierr);
      ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr);
      if (reason < 0  && reason != SNES_DIVERGED_MAX_IT) {
        snes->reason = SNES_DIVERGED_INNER;
        PetscFunctionReturn(0);
      }
      ierr = SNESGetFunction(snes->pc, &FPC, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr);
      ierr = VecCopy(FPC, F);CHKERRQ(ierr);
      ierr = SNESGetFunctionNorm(snes->pc, &fnorm);CHKERRQ(ierr);
    }

    /* Solve J Y = F, where J is Jacobian matrix */
    ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr);
    ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr);
    ierr = SNES_KSPSolve(snes,snes->ksp,F,Y);CHKERRQ(ierr);
    ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr);
    if (kspreason < 0) {
      if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
        ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr);
        snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
        break;
      }
    }
    ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr);
    snes->linear_its += lits;
    ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr);

    if (PetscLogPrintInfo){
      ierr = SNESNEWTONLSCheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr);
//.........这里部分代码省略.........

PetscErrorCode KSPSolve_NASH(KSP ksp)
{
#if defined(PETSC_USE_COMPLEX)
    SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "NASH is not available for complex systems");
#else
    KSP_NASH       *cg = (KSP_NASH*)ksp->data;
    PetscErrorCode ierr;
    MatStructure   pflag;
    Mat            Qmat, Mmat;
    Vec            r, z, p, d;
    PC             pc;

    PetscReal norm_r, norm_d, norm_dp1, norm_p, dMp;
    PetscReal alpha, beta, kappa, rz, rzm1;
    PetscReal rr, r2, step;

    PetscInt max_cg_its;

    PetscBool diagonalscale;

    PetscFunctionBegin;
    /***************************************************************************/
    /* Check the arguments and parameters.                                     */
    /***************************************************************************/

    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 (cg->radius < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");

    /***************************************************************************/
    /* Get the workspace vectors and initialize variables                      */
    /***************************************************************************/

    r2 = cg->radius * cg->radius;
    r  = ksp->work[0];
    z  = ksp->work[1];
    p  = ksp->work[2];
    d  = ksp->vec_sol;
    pc = ksp->pc;

    ierr = PCGetOperators(pc, &Qmat, &Mmat, &pflag);
    CHKERRQ(ierr);

    ierr       = VecGetSize(d, &max_cg_its);
    CHKERRQ(ierr);
    max_cg_its = PetscMin(max_cg_its, ksp->max_it);
    ksp->its   = 0;

    /***************************************************************************/
    /* Initialize objective function and direction.                            */
    /***************************************************************************/

    cg->o_fcn = 0.0;

    ierr       = VecSet(d, 0.0);
    CHKERRQ(ierr);            /* d = 0             */
    cg->norm_d = 0.0;

    /***************************************************************************/
    /* Begin the conjugate gradient method.  Check the right-hand side for     */
    /* numerical problems.  The check for not-a-number and infinite values     */
    /* need be performed only once.                                            */
    /***************************************************************************/

    ierr = VecCopy(ksp->vec_rhs, r);
    CHKERRQ(ierr);        /* r = -grad         */
    ierr = VecDot(r, r, &rr);
    CHKERRQ(ierr);               /* rr = r^T r        */
    if (PetscIsInfOrNanScalar(rr)) {
        /*************************************************************************/
        /* The right-hand side contains not-a-number or an infinite value.       */
        /* The gradient step does not work; return a zero value for the step.    */
        /*************************************************************************/

        ksp->reason = KSP_DIVERGED_NANORINF;
        ierr        = PetscInfo1(ksp, "KSPSolve_NASH: bad right-hand side: rr=%g\n", rr);
        CHKERRQ(ierr);
        PetscFunctionReturn(0);
    }

    /***************************************************************************/
    /* Check the preconditioner for numerical problems and for positive        */
    /* definiteness.  The check for not-a-number and infinite values need be   */
    /* performed only once.                                                    */
    /***************************************************************************/

    ierr = KSP_PCApply(ksp, r, z);
    CHKERRQ(ierr);          /* z = inv(M) r      */
    ierr = VecDot(r, z, &rz);
    CHKERRQ(ierr);               /* rz = r^T inv(M) r */
    if (PetscIsInfOrNanScalar(rz)) {
        /*************************************************************************/
        /* The preconditioner contains not-a-number or an infinite value.        */
        /* Return the gradient direction intersected with the trust region.      */
        /*************************************************************************/

        ksp->reason = KSP_DIVERGED_NANORINF;
        ierr        = PetscInfo1(ksp, "KSPSolve_NASH: bad preconditioner: rz=%g\n", rz);
        CHKERRQ(ierr);
//.........这里部分代码省略.........

PetscErrorCode SNESQNApply_LBFGS(SNES snes,PetscInt it,Vec Y,Vec X,Vec Xold,Vec D,Vec Dold)
{
  PetscErrorCode ierr;
  SNES_QN        *qn    = (SNES_QN*)snes->data;
  Vec            W      = snes->work[3];
  Vec            *dX    = qn->U;
  Vec            *dF    = qn->V;
  PetscScalar    *alpha = qn->alpha;
  PetscScalar    *beta  = qn->beta;
  PetscScalar    *dXtdF = qn->dXtdF;
  PetscScalar    *dFtdX = qn->dFtdX;
  PetscScalar    *YtdX  = qn->YtdX;

  /* ksp thing for Jacobian scaling */
  PetscInt           k,i,j,g,lits;
  PetscInt           m = qn->m;
  PetscScalar        t;
  PetscInt           l = m;

  PetscFunctionBegin;
  if (it < m) l = it;
  ierr = VecCopy(D,Y);CHKERRQ(ierr);
  if (it > 0) {
    k    = (it - 1) % l;
    ierr = VecCopy(D,dF[k]);CHKERRQ(ierr);
    ierr = VecAXPY(dF[k], -1.0, Dold);CHKERRQ(ierr);
    ierr = VecCopy(X, dX[k]);CHKERRQ(ierr);
    ierr = VecAXPY(dX[k], -1.0, Xold);CHKERRQ(ierr);
    if (qn->singlereduction) {
      ierr = VecMDotBegin(dF[k],l,dX,dXtdF);CHKERRQ(ierr);
      ierr = VecMDotBegin(dX[k],l,dF,dFtdX);CHKERRQ(ierr);
      ierr = VecMDotBegin(Y,l,dX,YtdX);CHKERRQ(ierr);
      ierr = VecMDotEnd(dF[k],l,dX,dXtdF);CHKERRQ(ierr);
      ierr = VecMDotEnd(dX[k],l,dF,dFtdX);CHKERRQ(ierr);
      ierr = VecMDotEnd(Y,l,dX,YtdX);CHKERRQ(ierr);
      for (j = 0; j < l; j++) {
        H(k, j) = dFtdX[j];
        H(j, k) = dXtdF[j];
      }
      /* copy back over to make the computation of alpha and beta easier */
      for (j = 0; j < l; j++) dXtdF[j] = H(j, j);
    } else {
      ierr = VecDot(dX[k], dF[k], &dXtdF[k]);CHKERRQ(ierr);
    }
    if (qn->scale_type == SNES_QN_SCALE_LINESEARCH) {
      ierr = SNESLineSearchGetLambda(snes->linesearch,&qn->scaling);CHKERRQ(ierr);
    }
  }

  /* outward recursion starting at iteration k's update and working back */
  for (i=0; i<l; i++) {
    k = (it-i-1)%l;
    if (qn->singlereduction) {
      /* construct t = dX[k] dot Y as Y_0 dot dX[k] + sum(-alpha[j]dX[k]dF[j]) */
      t = YtdX[k];
      for (j=0; j<i; j++) {
        g  = (it-j-1)%l;
        t -= alpha[g]*H(k, g);
      }
      alpha[k] = t/H(k,k);
    } else {
      ierr     = VecDot(dX[k],Y,&t);CHKERRQ(ierr);
      alpha[k] = t/dXtdF[k];
    }
    if (qn->monitor) {
      ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
      ierr = PetscViewerASCIIPrintf(qn->monitor, "it: %d k: %d alpha:        %14.12e\n", it, k, PetscRealPart(alpha[k]));CHKERRQ(ierr);
      ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
    }
    ierr = VecAXPY(Y,-alpha[k],dF[k]);CHKERRQ(ierr);
  }

  if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
    ierr = KSPSolve(snes->ksp,Y,W);CHKERRQ(ierr);
    SNESCheckKSPSolve(snes);
    ierr              = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr);
    snes->linear_its += lits;
    ierr              = VecCopy(W, Y);CHKERRQ(ierr);
  } else {
    ierr = VecScale(Y, qn->scaling);CHKERRQ(ierr);
  }
  if (qn->singlereduction) {
    ierr = VecMDot(Y,l,dF,YtdX);CHKERRQ(ierr);
  }
  /* inward recursion starting at the first update and working forward */
  for (i = 0; i < l; i++) {
    k = (it + i - l) % l;
    if (qn->singlereduction) {
      t = YtdX[k];
      for (j = 0; j < i; j++) {
        g  = (it + j - l) % l;
        t += (alpha[g] - beta[g])*H(g, k);
      }
      beta[k] = t / H(k, k);
    } else {
      ierr    = VecDot(dF[k], Y, &t);CHKERRQ(ierr);
      beta[k] = t / dXtdF[k];
    }
    ierr = VecAXPY(Y, (alpha[k] - beta[k]), dX[k]);CHKERRQ(ierr);
    if (qn->monitor) {
//.........这里部分代码省略.........

int main(int argc, char **argv)
{
  MPI_Comm          comm;
  PetscMPIInt       rank;
  PetscErrorCode    ierr;
  User              user;
  PetscLogDouble       v1, v2;
  PetscInt          nplot = 0;
  char              fileName[2048];


  ierr = PetscInitialize(&argc, &argv, (char*) 0, help);CHKERRQ(ierr);
  comm = PETSC_COMM_WORLD;
  ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
  ierr = PetscNew(&user);CHKERRQ(ierr);
  ierr = PetscNew(&user->algebra);CHKERRQ(ierr);
  ierr = PetscNew(&user->model);CHKERRQ(ierr);
  ierr = PetscNew(&user->model->physics);CHKERRQ(ierr);

  Algebra   algebra = user->algebra;

  ierr = LoadOptions(comm, user);CHKERRQ(ierr);
  ierr = PetscTime(&v1);CHKERRQ(ierr);
  ierr = CreateMesh(comm, user);CHKERRQ(ierr);
  ierr = PetscTime(&v2);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,
		       "Read and Distribute mesh takes %f sec \n", v2 - v1);CHKERRQ(ierr);
  ierr = SetUpLocalSpace(user);CHKERRQ(ierr); //Set up the dofs of each element
  ierr = ConstructGeometryFVM(&user->facegeom, &user->cellgeom, user);CHKERRQ(ierr);

  ierr = LimiterSetup(user);CHKERRQ(ierr);

  if (user->TimeIntegralMethod == EXPLICITMETHOD) { // explicit method
    if(user->myownexplicitmethod){// Using the fully explicit method based on my own routing
      ierr = PetscPrintf(PETSC_COMM_WORLD,"Using the fully explicit method based on my own routing\n");CHKERRQ(ierr);
      user->current_time = user->initial_time;
      user->current_step = 1;
      ierr = DMCreateGlobalVector(user->dm, &algebra->solution);CHKERRQ(ierr);
      ierr = PetscObjectSetName((PetscObject) algebra->solution, "solution");CHKERRQ(ierr);
      ierr = VecSet(algebra->solution, 0);CHKERRQ(ierr);
      ierr = SetInitialCondition(user->dm, algebra->solution, user);CHKERRQ(ierr);
      if(1){
        PetscViewer    viewer;
        ierr = OutputVTK(user->dm, "intialcondition.vtk", &viewer);CHKERRQ(ierr);
        ierr = VecView(algebra->solution, viewer);CHKERRQ(ierr);
        ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_WORLD,"Outputing the initial condition intialcondition.vtk!!! \n");CHKERRQ(ierr);
      }
      ierr = VecDuplicate(algebra->solution, &algebra->fn);CHKERRQ(ierr);
      ierr = VecDuplicate(algebra->solution, &algebra->oldsolution);CHKERRQ(ierr);
      if(user->Explicit_RK2||user->Explicit_RK4){
        ierr = PetscPrintf(PETSC_COMM_WORLD,"Use the second order Runge Kutta method \n");CHKERRQ(ierr);
      }else{
        ierr = PetscPrintf(PETSC_COMM_WORLD,"Use the first order forward Euler method \n");CHKERRQ(ierr);
      }
      nplot = 0; //the plot step
      while(user->current_time < (user->final_time - 0.05 * user->dt)){

        user->current_time = user->current_time + user->dt;
        ierr = FormTimeStepFunction(user, algebra, algebra->solution, algebra->fn);CHKERRQ(ierr);
        if(0){
            PetscViewer    viewer;
            ierr = OutputVTK(user->dm, "function.vtk", &viewer);CHKERRQ(ierr);
            ierr = VecView(algebra->fn, viewer);CHKERRQ(ierr);
            ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
          }

        if(user->Explicit_RK2){
        /*
          U^n_1   = U^n + 0.5*dt*f(U^n)
          U^{n+1} = U^n + dt*f(U^n_1)
        */
          ierr = VecCopy(algebra->solution, algebra->oldsolution);CHKERRQ(ierr);
          //note that algebra->oldsolution and algebra->solution are both U^n
          ierr = VecAXPY(algebra->solution, 0.5*user->dt, algebra->fn);CHKERRQ(ierr);
          //U^n_1 = U^n + 0.5*dt*f(U^n), now algebra->solution is U^n_1, and algebra->fn is f(U^n)

          ierr = FormTimeStepFunction(user, algebra, algebra->solution, algebra->fn);CHKERRQ(ierr);
          //algebra->fn is f(U^n_1)

          // reset the algebra->solution to U^n
          ierr = VecCopy(algebra->oldsolution, algebra->solution);CHKERRQ(ierr);
          ierr = VecAXPY(algebra->solution, user->dt, algebra->fn);CHKERRQ(ierr);
          // now algebra->solution is U^{n+1} = U^n + dt*f(U^n_1)
        }else if(user->Explicit_RK4){
        /* refer to https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods
          k_1     = f(U^n)
          U^n_1   = U^n + 0.5*dt*k_1
          k_2     = f(U^n_1)
          U^n_2   = U^n + 0.5*dt*k_2
          k_3     = f(U^n_2)
          U^n_3   = U^n + 0.5*dt*k_3
          k_4     = f(U^n_3)

          U^{n+1} = U^n + dt/6*(k_1 + 2*k_2 + 2*k_3 + k_4)
        */
          Vec  VecTemp; // store the U^n_1
          Vec  k1, k2, k3, k4;

          ierr = VecDuplicate(algebra->solution, &k1);CHKERRQ(ierr);
//.........这里部分代码省略.........

void CellBasedPdeHandler<DIM>::SolvePdeAndWriteResultsToFile(unsigned samplingTimestepMultiple)
{
    // Record whether we are solving PDEs on a coarse mesh
    bool using_coarse_pde_mesh = (mpCoarsePdeMesh != NULL);

    // If solving PDEs on a coarse mesh, each PDE should have an averaged source term; otherwise none should
    assert(!mPdeAndBcCollection.empty());
    for (unsigned pde_index=0; pde_index<mPdeAndBcCollection.size(); pde_index++)
    {
        assert(mPdeAndBcCollection[pde_index]);
        assert(mPdeAndBcCollection[pde_index]->HasAveragedSourcePde() == using_coarse_pde_mesh || dynamic_cast<MultipleCaBasedCellPopulation<DIM>*>(mpCellPopulation));
    }

    // Make sure the cell population is in a nice state
    mpCellPopulation->Update();

    // Store a pointer to the (population-level or coarse) mesh
    TetrahedralMesh<DIM,DIM>* p_mesh;
    if (using_coarse_pde_mesh)
    {
        p_mesh = mpCoarsePdeMesh;
    }
    else
    {
        // If not using a coarse PDE mesh, we must be using a MeshBasedCellPopulation
        p_mesh = &(static_cast<MeshBasedCellPopulation<DIM>*>(mpCellPopulation)->rGetMesh());
    }

    // Loop over elements of mPdeAndBcCollection
    for (unsigned pde_index=0; pde_index<mPdeAndBcCollection.size(); pde_index++)
    {
        // Get pointer to this PdeAndBoundaryConditions object
        PdeAndBoundaryConditions<DIM>* p_pde_and_bc = mPdeAndBcCollection[pde_index];

        // Set up boundary conditions
        std::auto_ptr<BoundaryConditionsContainer<DIM,DIM,1> > p_bcc = ConstructBoundaryConditionsContainer(p_pde_and_bc, p_mesh);

        // If the solution at the previous timestep exists...
        PetscInt previous_solution_size = 0;
        if (p_pde_and_bc->GetSolution())
        {
            VecGetSize(p_pde_and_bc->GetSolution(), &previous_solution_size);
        }

        // ...then record whether it is the correct size...
        bool is_previous_solution_size_correct = (previous_solution_size == (int)p_mesh->GetNumNodes());

        // ...and if it is, store it as an initial guess for the PDE solver
        Vec initial_guess;
        if (is_previous_solution_size_correct)
        {
            // This Vec is copied by the solver's Solve() method, so must be deleted here too
            VecDuplicate(p_pde_and_bc->GetSolution(), &initial_guess);
            VecCopy(p_pde_and_bc->GetSolution(), initial_guess);
            p_pde_and_bc->DestroySolution();
        }
        else
        {
            ///\todo enable the coarse PDE mesh to change size, e.g. for a growing domain (#630/#1891)
            if (!using_coarse_pde_mesh && p_pde_and_bc->GetSolution())
            {
                assert(previous_solution_size != 0);
                p_pde_and_bc->DestroySolution();
            }
        }

        // Create a PDE solver and solve the PDE on the (population-level or coarse) mesh
        if (p_pde_and_bc->HasAveragedSourcePde())
        {
            // When using a coarse PDE mesh, we must set up the source terms before solving the PDE.
            // Pass in mCellPdeElementMap to speed up finding cells.
            this->UpdateCellPdeElementMap();
            p_pde_and_bc->SetUpSourceTermsForAveragedSourcePde(p_mesh, &mCellPdeElementMap);


            SimpleLinearEllipticSolver<DIM,DIM> solver(p_mesh, p_pde_and_bc->GetPde(), p_bcc.get());

            // If we have an initial guess, use this when solving the system...
            if (is_previous_solution_size_correct)
            {
                p_pde_and_bc->SetSolution(solver.Solve(initial_guess));
                PetscTools::Destroy(initial_guess);
            }
            else // ...otherwise do not supply one
            {
                p_pde_and_bc->SetSolution(solver.Solve());
            }
        }
        else
        {
            CellBasedPdeSolver<DIM> solver(p_mesh, p_pde_and_bc->GetPde(), p_bcc.get());

            // If we have an initial guess, use this...
            if (is_previous_solution_size_correct)
            {
                p_pde_and_bc->SetSolution(solver.Solve(initial_guess));
                PetscTools::Destroy(initial_guess);
            }
            else // ...otherwise do not supply one
            {
//.........这里部分代码省略.........

int main(int argc,char **argv)
{
  Vec            x1, x2, y1, y2, y3, y4;
  PetscViewer    viewer;
  PetscMPIInt    rank;
  PetscInt       i, nlocal, n = 6;
  PetscScalar   *array;
  PetscBool      equal;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = PetscInitialize(&argc, &argv, (char *) 0, help);CHKERRQ(ierr); 
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRQ(ierr);
  ierr = PetscOptionsGetInt(PETSC_NULL, "-n", &n, PETSC_NULL);CHKERRQ(ierr);

  ierr = VecCreate(PETSC_COMM_WORLD, &x1);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) x1, "TestVec");CHKERRQ(ierr);
  ierr = VecSetSizes(x1, PETSC_DECIDE, n);CHKERRQ(ierr);
  ierr = VecSetFromOptions(x1);CHKERRQ(ierr);

  ierr = VecGetLocalSize(x1, &nlocal);CHKERRQ(ierr);
  ierr = VecGetArray(x1, &array);CHKERRQ(ierr);
  for(i = 0; i < nlocal; i++) {
    array[i] = rank + 1;
  }
  ierr = VecRestoreArray(x1, &array);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(x1);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(x1);CHKERRQ(ierr);

  ierr = VecCreate(PETSC_COMM_WORLD, &x2);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) x2, "TestVec2");CHKERRQ(ierr);
  ierr = VecSetSizes(x2, PETSC_DECIDE, n);CHKERRQ(ierr);
  ierr = VecSetBlockSize(x2, 2);CHKERRQ(ierr);
  ierr = VecSetFromOptions(x2);CHKERRQ(ierr);
  ierr = VecCopy(x1, x2);CHKERRQ(ierr);

  ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD, "ex19.h5", FILE_MODE_WRITE, &viewer);CHKERRQ(ierr);
  ierr = PetscViewerHDF5PushGroup(viewer, "/");CHKERRQ(ierr);
  ierr = VecView(x1, viewer);CHKERRQ(ierr);
  ierr = PetscViewerHDF5PushGroup(viewer, "/testBlockSize");CHKERRQ(ierr);
  ierr = VecView(x2, viewer);CHKERRQ(ierr);
  ierr = PetscViewerHDF5PushGroup(viewer, "/testTimestep");CHKERRQ(ierr);
  ierr = PetscViewerHDF5SetTimestep(viewer, 0);CHKERRQ(ierr);
  ierr = VecView(x2, viewer);CHKERRQ(ierr);
  ierr = PetscViewerHDF5SetTimestep(viewer, 1);CHKERRQ(ierr);
  ierr = VecView(x2, viewer);CHKERRQ(ierr);
  ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
  ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
  ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
  ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);

  ierr = VecCreate(PETSC_COMM_WORLD, &y1);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) y1, "TestVec");CHKERRQ(ierr);
  ierr = VecSetSizes(y1, PETSC_DECIDE, n);CHKERRQ(ierr);
  ierr = VecSetFromOptions(y1);CHKERRQ(ierr);

  ierr = VecCreate(PETSC_COMM_WORLD,&y2);CHKERRQ(ierr);
  ierr = VecSetBlockSize(y2, 2);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) y2, "TestVec2");CHKERRQ(ierr);

  ierr = VecCreate(PETSC_COMM_WORLD,&y3);CHKERRQ(ierr);
  ierr = VecSetBlockSize(y3, 2);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) y3, "TestVec2");CHKERRQ(ierr);

  ierr = VecCreate(PETSC_COMM_WORLD,&y4);CHKERRQ(ierr);
  ierr = VecSetBlockSize(y4, 2);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) y4, "TestVec2");CHKERRQ(ierr);

  ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD, "ex19.h5", FILE_MODE_READ, &viewer);CHKERRQ(ierr);
  ierr = PetscViewerHDF5PushGroup(viewer, "/");CHKERRQ(ierr);
  ierr = VecLoad(y1, viewer);CHKERRQ(ierr);
  
  ierr = PetscViewerHDF5PushGroup(viewer, "/testBlockSize");CHKERRQ(ierr);
  ierr = VecLoad(y2, viewer);CHKERRQ(ierr);
  ierr = PetscViewerHDF5PushGroup(viewer, "/testTimestep");CHKERRQ(ierr);
  ierr = PetscViewerHDF5SetTimestep(viewer, 0);CHKERRQ(ierr);
  ierr = VecLoad(y3, viewer);CHKERRQ(ierr);
  ierr = PetscViewerHDF5SetTimestep(viewer, 1);CHKERRQ(ierr);
  ierr = VecLoad(y4, viewer);CHKERRQ(ierr);
  ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
  ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
  ierr = PetscViewerHDF5PopGroup(viewer);CHKERRQ(ierr);
  ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);

  ierr = VecEqual(x1, y1, &equal);CHKERRQ(ierr);
  if (!equal) {
    ierr = VecView(x1, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
    ierr = VecView(y1, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
    SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Error in HDF5 viewer");
  }
  ierr = VecEqual(x2, y2, &equal);CHKERRQ(ierr);
  if (!equal) {
    ierr = VecView(x2, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
    ierr = VecView(y2, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
    SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Error in HDF5 viewer");
  }

  ierr = VecDestroy(&x1);CHKERRQ(ierr);
  ierr = VecDestroy(&x2);CHKERRQ(ierr);
  ierr = VecDestroy(&y1);CHKERRQ(ierr);
//.........这里部分代码省略.........

PetscErrorCode  KSPSolve_GROPPCG(KSP ksp)
{
  PetscErrorCode ierr;
  PetscInt       i;
  PetscScalar    alpha,beta = 0.0,gamma,gammaNew,t;
  PetscReal      dp = 0.0;
  Vec            x,b,r,p,s,S,z,Z;
  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);

  x = ksp->vec_sol;
  b = ksp->vec_rhs;
  r = ksp->work[0];
  p = ksp->work[1];
  s = ksp->work[2];
  S = ksp->work[3];
  z = ksp->work[4];
  Z = ksp->work[5];

  ierr = PCGetOperators(ksp->pc,&Amat,&Pmat);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) */
  }

  ierr = KSP_PCApply(ksp,r,z);CHKERRQ(ierr);                   /*     z <- Br   */
  ierr = VecCopy(z,p);CHKERRQ(ierr);                           /*     p <- z    */
  ierr = VecDotBegin(r,z,&gamma);CHKERRQ(ierr);                  /*     gamma <- z'*r       */
  ierr = PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)r));CHKERRQ(ierr);
  ierr = KSP_MatMult(ksp,Amat,p,s);CHKERRQ(ierr);              /*     s <- Ap   */
  ierr = VecDotEnd(r,z,&gamma);CHKERRQ(ierr);                  /*     gamma <- z'*r       */

  switch (ksp->normtype) {
  case KSP_NORM_PRECONDITIONED:
    /* This could be merged with the computation of gamma above */
    ierr = VecNorm(z,NORM_2,&dp);CHKERRQ(ierr);                /*     dp <- z'*z = e'*A'*B'*B*A'*e'     */
    break;
  case KSP_NORM_UNPRECONDITIONED:
    /* This could be merged with the computation of gamma above */
    ierr = VecNorm(r,NORM_2,&dp);CHKERRQ(ierr);                /*     dp <- r'*r = e'*A'*A*e            */
    break;
  case KSP_NORM_NATURAL:
    if (PetscIsInfOrNanScalar(gamma)) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_FP,"Infinite or not-a-number generated in dot product");
    dp = PetscSqrtReal(PetscAbsScalar(gamma));                  /*     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);

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

    ierr = VecDotBegin(p,s,&t);CHKERRQ(ierr);
    ierr = PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)p));CHKERRQ(ierr);

    ierr = KSP_PCApply(ksp,s,S);CHKERRQ(ierr);         /*   S <- Bs       */

    ierr = VecDotEnd(p,s,&t);CHKERRQ(ierr);

    alpha = gamma / t;
    ierr  = VecAXPY(x, alpha,p);CHKERRQ(ierr);   /*     x <- x + alpha * p   */
    ierr  = VecAXPY(r,-alpha,s);CHKERRQ(ierr);   /*     r <- r - alpha * s   */
    ierr  = VecAXPY(z,-alpha,S);CHKERRQ(ierr);   /*     z <- z - alpha * S   */

    if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
      ierr = VecNormBegin(r,NORM_2,&dp);CHKERRQ(ierr);
    } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
      ierr = VecNormBegin(z,NORM_2,&dp);CHKERRQ(ierr);
    }
    ierr = VecDotBegin(r,z,&gammaNew);CHKERRQ(ierr);
    ierr = PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)r));CHKERRQ(ierr);

    ierr = KSP_MatMult(ksp,Amat,z,Z);CHKERRQ(ierr);      /*   Z <- Az       */

    if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
      ierr = VecNormEnd(r,NORM_2,&dp);CHKERRQ(ierr);
    } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
      ierr = VecNormEnd(z,NORM_2,&dp);CHKERRQ(ierr);
    }
    ierr = VecDotEnd(r,z,&gammaNew);CHKERRQ(ierr);

    if (ksp->normtype == KSP_NORM_NATURAL) {
      if (PetscIsInfOrNanScalar(gammaNew)) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_FP,"Infinite or not-a-number generated in dot product");
//.........这里部分代码省略.........

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 MatMultAdd_SeqSBSTRM_5(Mat A,Vec xx,Vec yy,Vec zz)
{
  Mat_SeqSBAIJ   *a      = (Mat_SeqSBAIJ*)A->data;
  Mat_SeqSBSTRM  *sbstrm = (Mat_SeqSBSTRM*)A->spptr;
  PetscScalar    *x,*xb, *z;
  MatScalar      *v1, *v2, *v3, *v4, *v5;
  PetscScalar    x1, x2, x3, x4, x5;
  PetscScalar    xr1, xr2, xr3, xr4, xr5;
  PetscScalar    sum1, sum2, sum3, sum4, sum5;
  PetscErrorCode ierr;
  PetscInt       mbs       =a->mbs,i,*aj=a->j,*ai=a->i,n,*ib,cval,j,jmin;
  PetscInt       nonzerorow=0;
  PetscInt       slen;

  PetscFunctionBegin;
  ierr = VecCopy(yy,zz);CHKERRQ(ierr);
  ierr = VecGetArray(xx,&x);CHKERRQ(ierr);
  ierr = VecGetArray(zz,&z);CHKERRQ(ierr);


  slen = 5*(ai[mbs]-ai[0]);
  v1   = sbstrm->as;
  v2   = v1 + slen;
  v3   = v2 + slen;
  v4   = v3 + slen;
  v5   = v4 + slen;

  xb = x;

  for (i=0; i<mbs; i++) {
    n  = ai[i+1] - ai[i];
    x1 = xb[0]; x2 = xb[1]; x3 = xb[2]; x4 = xb[3]; x5 = xb[4]; xb += 5;

    sum1 = z[5*i];
    sum2 = z[5*i+1];
    sum3 = z[5*i+2];
    sum4 = z[5*i+3];
    sum5 = z[5*i+4];

    nonzerorow += (n>0);

    jmin = 0;
    ib   = aj + ai[i];
    if (*ib == i) {     /* (diag of A)*x, only upper triangular part is used */
      sum1 += v1[0]*x1 + v1[1]*x2 + v1[2]*x3 + v1[3]*x4 + v1[4]*x5;
      sum2 += v1[1]*x1 + v2[1]*x2 + v2[2]*x3 + v2[3]*x4 + v2[4]*x5;
      sum3 += v1[2]*x1 + v2[2]*x2 + v3[2]*x3 + v3[3]*x4 + v3[4]*x5;
      sum4 += v1[3]*x1 + v2[3]*x2 + v3[3]*x3 + v4[3]*x4 + v4[4]*x5;
      sum5 += v1[4]*x1 + v2[4]*x2 + v3[4]*x3 + v4[4]*x4 + v5[4]*x5;

      v1 += 5; v2 += 5; v3 += 5; v4 += 5; v5 += 5;
      jmin++;
    }

    for (j=jmin; j<n; j++) {
      cval       = ib[j]*5;
      z[cval]   += v1[0]*x1 + v2[0]*x2 + v3[0]*x3 + v4[0]*x4 + v5[0]*x5;
      z[cval+1] += v1[1]*x1 + v2[1]*x2 + v3[1]*x3 + v4[1]*x4 + v5[1]*x5;
      z[cval+2] += v1[2]*x1 + v2[2]*x2 + v3[2]*x3 + v4[2]*x4 + v5[2]*x5;
      z[cval+3] += v1[3]*x1 + v2[3]*x2 + v3[3]*x3 + v4[3]*x4 + v5[3]*x5;
      z[cval+4] += v1[4]*x1 + v2[4]*x2 + v3[4]*x3 + v4[4]*x4 + v5[4]*x5;

      xr1 = x[cval]; xr2 = x[cval+1]; xr3 = x[cval+2]; xr4 = x[cval+3]; xr5 = x[cval+4];

      sum1 += v1[0]*xr1 + v1[1]*xr2 + v1[2]*xr3 + v1[3]*xr4 + v1[4]*xr5;
      sum2 += v2[0]*xr1 + v2[1]*xr2 + v2[2]*xr3 + v2[3]*xr4 + v2[4]*xr5;
      sum3 += v3[0]*xr1 + v3[1]*xr2 + v3[2]*xr3 + v3[3]*xr4 + v3[4]*xr5;
      sum4 += v4[0]*xr1 + v4[1]*xr2 + v4[2]*xr3 + v4[3]*xr4 + v4[4]*xr5;
      sum5 += v5[0]*xr1 + v5[1]*xr2 + v5[2]*xr3 + v5[3]*xr4 + v5[4]*xr5;

      v1 += 5; v2 += 5; v3 += 5; v4 += 5; v5 += 5;
    }
    z[5*i]   = sum1;
    z[5*i+1] = sum2;
    z[5*i+2] = sum3;
    z[5*i+3] = sum4;
    z[5*i+4] = sum5;
  }
  ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr);
  ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr);
  ierr = PetscLogFlops(50.0*(a->nz*2.0 - nonzerorow));CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

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 */
  loc_it = 0;

  /* LGMRES_MOD: determine number of arnoldi steps to take */
  /* if approx_constant then we keep the space the same size even if
     we don't have the full number of aug vectors yet*/
  if (lgmres->approx_constant) it_arnoldi = max_k - lgmres->aug_ct;
  else it_arnoldi = max_k - aug_dim;

  it_total =  it_arnoldi + lgmres->aug_ct;

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

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

  /* check for the convergence */
  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);
  }

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

  ksp->rnorm = res;


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


  /* 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 */
  ierr = (*ksp->converged)(ksp,ksp->its,res,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);

  while (!ksp->reason && loc_it < it_total && ksp->its < max_it) { /* LGMRES_MOD: changed to it_total */
    ierr       = KSPLogResidualHistory(ksp,res);CHKERRQ(ierr);
    lgmres->it = (loc_it - 1);
    ierr       = KSPMonitor(ksp,ksp->its,res);CHKERRQ(ierr);

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

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

static PetscErrorCode SNESLineSearchApply_CP(SNESLineSearch linesearch)
{
  PetscBool      changed_y, changed_w;
  PetscErrorCode ierr;
  Vec            X, Y, F, W;
  SNES           snes;
  PetscReal      xnorm, ynorm, gnorm, steptol, atol, rtol, ltol, maxstep;

  PetscReal   lambda, lambda_old, lambda_update, delLambda;
  PetscScalar fty, fty_init, fty_old, fty_mid1, fty_mid2, s;
  PetscInt    i, max_its;

  PetscViewer monitor;

  PetscFunctionBegin;
  ierr = SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, NULL);CHKERRQ(ierr);
  ierr = SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm);CHKERRQ(ierr);
  ierr = SNESLineSearchGetSNES(linesearch, &snes);CHKERRQ(ierr);
  ierr = SNESLineSearchGetLambda(linesearch, &lambda);CHKERRQ(ierr);
  ierr = SNESLineSearchGetTolerances(linesearch, &steptol, &maxstep, &rtol, &atol, &ltol, &max_its);CHKERRQ(ierr);
  ierr = SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED);CHKERRQ(ierr);
  ierr = SNESLineSearchGetDefaultMonitor(linesearch, &monitor);CHKERRQ(ierr);

  /* precheck */
  ierr       = SNESLineSearchPreCheck(linesearch,X,Y,&changed_y);CHKERRQ(ierr);
  lambda_old = 0.0;

  ierr = VecDot(F,Y,&fty_old);CHKERRQ(ierr);
  fty_init = fty_old;

  for (i = 0; i < max_its; i++) {
    /* compute the norm at lambda */
    ierr = VecCopy(X, W);CHKERRQ(ierr);
    ierr = VecAXPY(W, -lambda, Y);CHKERRQ(ierr);
    if (linesearch->ops->viproject) {
      ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
    }
    ierr = (*linesearch->ops->snesfunc)(snes,W,F);CHKERRQ(ierr);
    ierr = VecDot(F,Y,&fty);CHKERRQ(ierr);

    delLambda = lambda - lambda_old;

    /* check for convergence */
    if (PetscAbsReal(delLambda) < steptol*lambda) break;
    if (PetscAbsScalar(fty) / PetscAbsScalar(fty_init) < rtol) break;
    if (PetscAbsScalar(fty) < atol && i > 0) break;
    if (monitor) {
      ierr = PetscViewerASCIIAddTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
      ierr = PetscViewerASCIIPrintf(monitor,"    Line search: lambdas = [%g, %g], ftys = [%g, %g]\n",(double)lambda, (double)lambda_old, (double)PetscRealPart(fty), (double)PetscRealPart(fty_old));CHKERRQ(ierr);
      ierr = PetscViewerASCIISubtractTab(monitor,((PetscObject)linesearch)->tablevel);CHKERRQ(ierr);
    }

    /* compute the search direction */
    if (linesearch->order == SNES_LINESEARCH_ORDER_LINEAR) {
      s = (fty - fty_old) / delLambda;
    } else if (linesearch->order == SNES_LINESEARCH_ORDER_QUADRATIC) {
      ierr = VecCopy(X, W);CHKERRQ(ierr);
      ierr = VecAXPY(W, -0.5*(lambda + lambda_old), Y);CHKERRQ(ierr);
      if (linesearch->ops->viproject) {
        ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
      }
      ierr = (*linesearch->ops->snesfunc)(snes,W,F);CHKERRQ(ierr);
      ierr = VecDot(F, Y, &fty_mid1);CHKERRQ(ierr);
      s    = (3.*fty - 4.*fty_mid1 + fty_old) / delLambda;
    } else {
      ierr = VecCopy(X, W);CHKERRQ(ierr);
      ierr = VecAXPY(W, -0.5*(lambda + lambda_old), Y);CHKERRQ(ierr);
      if (linesearch->ops->viproject) {
        ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
      }
      ierr = (*linesearch->ops->snesfunc)(snes,W,F);CHKERRQ(ierr);
      ierr = VecDot(F, Y, &fty_mid1);CHKERRQ(ierr);
      ierr = VecCopy(X, W);CHKERRQ(ierr);
      ierr = VecAXPY(W, -(lambda + 0.5*(lambda - lambda_old)), Y);CHKERRQ(ierr);
      if (linesearch->ops->viproject) {
        ierr = (*linesearch->ops->viproject)(snes, W);CHKERRQ(ierr);
      }
      ierr = (*linesearch->ops->snesfunc)(snes, W, F);CHKERRQ(ierr);
      ierr = VecDot(F, Y, &fty_mid2);CHKERRQ(ierr);
      s    = (2.*fty_mid2 + 3.*fty - 6.*fty_mid1 + fty_old) / (3.*delLambda);
    }
    /* if the solve is going in the wrong direction, fix it */
    if (PetscRealPart(s) > 0.) s = -s;
    lambda_update =  lambda - PetscRealPart(fty / s);

    /* switch directions if we stepped out of bounds */
    if (lambda_update < steptol) lambda_update = lambda + PetscRealPart(fty / s);

    if (PetscIsInfOrNanReal(lambda_update)) break;
    if (lambda_update > maxstep) break;

    /* compute the new state of the line search */
    lambda_old = lambda;
    lambda     = lambda_update;
    fty_old    = fty;
  }
  /* construct the solution */
  ierr = VecCopy(X, W);CHKERRQ(ierr);
  ierr = VecAXPY(W, -lambda, Y);CHKERRQ(ierr);
  if (linesearch->ops->viproject) {
//.........这里部分代码省略.........

static PetscErrorCode TaoSolve_SQPCON(Tao tao)
{
  TAO_SQPCON                   *sqpconP = (TAO_SQPCON*)tao->data;
  PetscInt                     iter=0;
  TaoConvergedReason           reason = TAO_CONTINUE_ITERATING;
  TaoLineSearchConvergedReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING;
  PetscReal                    step=1.0,f,fm, fold;
  PetscReal                    cnorm, mnorm;
  PetscBool                    use_update=PETSC_TRUE; /*  don't update Q if line search failed */
  PetscErrorCode               ierr;

  PetscFunctionBegin;
  /* Scatter to U,V */
  ierr = VecScatterBegin(sqpconP->state_scatter, tao->solution, sqpconP->U, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr);
  ierr = VecScatterEnd(sqpconP->state_scatter, tao->solution, sqpconP->U, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr);
  ierr = VecScatterBegin(sqpconP->design_scatter, tao->solution, sqpconP->V, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr);
  ierr = VecScatterEnd(sqpconP->design_scatter, tao->solution, sqpconP->V, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr);

  /* Evaluate Function, Gradient, Constraints, and Jacobian */
  ierr = TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);CHKERRQ(ierr);
  ierr = TaoComputeConstraints(tao,tao->solution, tao->constraints);CHKERRQ(ierr);
  ierr = TaoComputeJacobianState(tao,tao->solution, &tao->jacobian_state, &tao->jacobian_state_pre, &tao->jacobian_state_inv, &sqpconP->statematflag);CHKERRQ(ierr);
  ierr = TaoComputeJacobianDesign(tao,tao->solution, &tao->jacobian_design, &tao->jacobian_design_pre, &sqpconP->statematflag);CHKERRQ(ierr);

  /* Scatter gradient to GU,GV */
  ierr = VecScatterBegin(sqpconP->state_scatter, tao->gradient, sqpconP->GU, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr);
  ierr = VecScatterEnd(sqpconP->state_scatter, tao->gradient, sqpconP->GU, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr);
  ierr = VecScatterBegin(sqpconP->design_scatter, tao->gradient, sqpconP->GV, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr);
  ierr = VecScatterEnd(sqpconP->design_scatter, tao->gradient, sqpconP->GV, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr);
  ierr = VecNorm(tao->gradient, NORM_2, &mnorm);CHKERRQ(ierr);

  /* Evaluate constraint norm */
  ierr = VecNorm(tao->constraints, NORM_2, &cnorm);CHKERRQ(ierr);

  /* Monitor convergence */
  ierr = TaoMonitor(tao, iter,f,mnorm,cnorm,step,&reason);CHKERRQ(ierr);

  while (reason == TAO_CONTINUE_ITERATING) {
    /* Solve tbar = -A\t (t is constraints vector) */
    ierr = MatMult(tao->jacobian_state_inv, tao->constraints, sqpconP->Tbar);CHKERRQ(ierr);
    ierr = VecScale(sqpconP->Tbar, -1.0);CHKERRQ(ierr);

    /* aqwac =  A'\(Q*Tbar + c) */
    if (iter > 0) {
      ierr = MatMult(sqpconP->Q,sqpconP->Tbar,sqpconP->WV);CHKERRQ(ierr);
    } else {
      ierr = VecCopy(sqpconP->Tbar, sqpconP->WV);CHKERRQ(ierr);
    }
    ierr = VecAXPY(sqpconP->WV,1.0,sqpconP->GU);CHKERRQ(ierr);

    ierr = MatMultTranspose(tao->jacobian_state_inv, sqpconP->WV, sqpconP->aqwac);CHKERRQ(ierr);

    /* Reduced Gradient dbar = d -  B^t * aqwac */
    ierr = MatMultTranspose(tao->jacobian_design,sqpconP->aqwac, sqpconP->dbar);CHKERRQ(ierr);
    ierr = VecScale(sqpconP->dbar, -1.0);CHKERRQ(ierr);
    ierr = VecAXPY(sqpconP->dbar,1.0,sqpconP->GV);CHKERRQ(ierr);

    /* update reduced hessian */
    ierr = MatLMVMUpdate(sqpconP->R, sqpconP->V, sqpconP->dbar);CHKERRQ(ierr);

    /* Solve R*dv = -dbar using approx. hessian */
    ierr = MatLMVMSolve(sqpconP->R, sqpconP->dbar, sqpconP->DV);CHKERRQ(ierr);
    ierr = VecScale(sqpconP->DV, -1.0);CHKERRQ(ierr);

    /* Backsolve for u =  A\(g - B*dv)  = tbar - A\(B*dv)*/
    ierr = MatMult(tao->jacobian_design, sqpconP->DV, sqpconP->WL);CHKERRQ(ierr);
    ierr = MatMult(tao->jacobian_state_inv, sqpconP->WL, sqpconP->DU);CHKERRQ(ierr);
    ierr = VecScale(sqpconP->DU, -1.0);CHKERRQ(ierr);
    ierr = VecAXPY(sqpconP->DU, 1.0, sqpconP->Tbar);CHKERRQ(ierr);

    /* Assemble Big D */
    ierr = VecScatterBegin(sqpconP->state_scatter, sqpconP->DU, tao->stepdirection, INSERT_VALUES, SCATTER_REVERSE);CHKERRQ(ierr);
    ierr = VecScatterEnd(sqpconP->state_scatter, sqpconP->DU, tao->stepdirection, INSERT_VALUES, SCATTER_REVERSE);CHKERRQ(ierr);
    ierr = VecScatterBegin(sqpconP->design_scatter, sqpconP->DV, tao->stepdirection, INSERT_VALUES, SCATTER_REVERSE);CHKERRQ(ierr);
    ierr = VecScatterEnd(sqpconP->design_scatter, sqpconP->DV, tao->stepdirection, INSERT_VALUES, SCATTER_REVERSE);CHKERRQ(ierr);

    /* Perform Line Search */
    ierr = VecCopy(tao->solution, sqpconP->Xold);CHKERRQ(ierr);
    ierr = VecCopy(tao->gradient, sqpconP->Gold);CHKERRQ(ierr);
    fold = f;
    ierr = TaoLineSearchComputeObjectiveAndGradient(tao->linesearch,tao->solution,&fm,sqpconP->GL);CHKERRQ(ierr);
    ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
    ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &fm, sqpconP->GL, tao->stepdirection,&step, &ls_reason);CHKERRQ(ierr);
    ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr);
    if (ls_reason < 0) {
      ierr = VecCopy(sqpconP->Xold, tao->solution);
      ierr = VecCopy(sqpconP->Gold, tao->gradient);
      f = fold;
      ierr = VecAXPY(tao->solution, 1.0, tao->stepdirection);CHKERRQ(ierr);
      ierr = PetscInfo(tao,"Line Search Failed, using full step.");CHKERRQ(ierr);
      use_update=PETSC_FALSE;
    } else {
      use_update = PETSC_TRUE;
    }

    /* Scatter X to U,V */
    ierr = VecScatterBegin(sqpconP->state_scatter, tao->solution, sqpconP->U, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr);
    ierr = VecScatterEnd(sqpconP->state_scatter, tao->solution, sqpconP->U, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr);
    ierr = VecScatterBegin(sqpconP->design_scatter, tao->solution, sqpconP->V, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr);
    ierr = VecScatterEnd(sqpconP->design_scatter, tao->solution, sqpconP->V, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      (*pc->ops->setfromoptions)(pc);
    }
    ierr = PCShellSetName(pc, "bfgs");CHKERRQ(ierr);
    ierr = PCShellSetContext(pc, tr->M);CHKERRQ(ierr);
    ierr = PCShellSetApply(pc, MatLMVMSolveShell);CHKERRQ(ierr);
    break;

  default:
    /*  Use the pc method set by pc_type */
    break;
  }

  /*  Initialize trust-region radius */
  switch(tr->init_type) {
  case NTR_INIT_CONSTANT:
    /*  Use the initial radius specified */
    break;

  case NTR_INIT_INTERPOLATION:
    /*  Use the initial radius specified */
    max_radius = 0.0;

    for (j = 0; j < j_max; ++j) {
      fmin = f;
      sigma = 0.0;

      if (needH) {
        ierr = TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);CHKERRQ(ierr);
        needH = 0;
      }

      for (i = 0; i < i_max; ++i) {

        ierr = VecCopy(tao->solution, tr->W);CHKERRQ(ierr);
        ierr = VecAXPY(tr->W, -tao->trust/gnorm, tao->gradient);CHKERRQ(ierr);
        ierr = TaoComputeObjective(tao, tr->W, &ftrial);CHKERRQ(ierr);

        if (PetscIsInfOrNanReal(ftrial)) {
          tau = tr->gamma1_i;
        }
        else {
          if (ftrial < fmin) {
            fmin = ftrial;
            sigma = -tao->trust / gnorm;
          }

          ierr = MatMult(tao->hessian, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
          ierr = VecDot(tao->gradient, tao->stepdirection, &prered);CHKERRQ(ierr);

          prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm));
          actred = f - ftrial;
          if ((PetscAbsScalar(actred) <= tr->epsilon) &&
              (PetscAbsScalar(prered) <= tr->epsilon)) {
            kappa = 1.0;
          }
          else {
            kappa = actred / prered;
          }

          tau_1 = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust + (1.0 - tr->theta_i) * prered - actred);
          tau_2 = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust - (1.0 + tr->theta_i) * prered + actred);
          tau_min = PetscMin(tau_1, tau_2);
          tau_max = PetscMax(tau_1, tau_2);

          if (PetscAbsScalar(kappa - 1.0) <= tr->mu1_i) {
            /*  Great agreement */