/*
   PostSetSubKSP - Optional user-defined routine that reset SubKSP options when hierarchical bjacobi PC is used
   e.g,
     mpiexec -n 8 ./ex3 -nox -n 10000 -ksp_type fgmres -pc_type bjacobi -pc_bjacobi_blocks 4 -sub_ksp_type gmres -sub_ksp_max_it 3 -post_setsubksp -sub_ksp_rtol 1.e-16
   Set by SNESLineSearchSetPostCheck().

   Input Parameters:
   linesearch - the LineSearch context
   xcurrent - current solution
   y - search direction and length
   x    - the new candidate iterate

   Output Parameters:
   y    - proposed step (search direction and length) (possibly changed)
   x    - current iterate (possibly modified)

 */
PetscErrorCode PostSetSubKSP(SNESLineSearch linesearch,Vec xcurrent,Vec y,Vec x,PetscBool  *changed_y,PetscBool  *changed_x, void * ctx)
{
  PetscErrorCode ierr;
  SetSubKSPCtx   *check;
  PetscInt       iter,its,sub_its,maxit;
  KSP            ksp,sub_ksp,*sub_ksps;
  PC             pc;
  PetscReal      ksp_ratio;
  SNES           snes;

  PetscFunctionBeginUser;
  ierr    = SNESLineSearchGetSNES(linesearch, &snes);CHKERRQ(ierr);
  check   = (SetSubKSPCtx*)ctx;
  ierr    = SNESGetIterationNumber(snes,&iter);CHKERRQ(ierr);
  ierr    = SNESGetKSP(snes,&ksp);CHKERRQ(ierr);
  ierr    = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
  ierr    = PCBJacobiGetSubKSP(pc,NULL,NULL,&sub_ksps);CHKERRQ(ierr);
  sub_ksp = sub_ksps[0];
  ierr    = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);      /* outer KSP iteration number */
  ierr    = KSPGetIterationNumber(sub_ksp,&sub_its);CHKERRQ(ierr); /* inner KSP iteration number */

  if (iter) {
    ierr      = PetscPrintf(PETSC_COMM_WORLD,"    ...PostCheck snes iteration %D, ksp_it %d %d, subksp_it %d\n",iter,check->its0,its,sub_its);CHKERRQ(ierr);
    ksp_ratio = ((PetscReal)(its))/check->its0;
    maxit     = (PetscInt)(ksp_ratio*sub_its + 0.5);
    if (maxit < 2) maxit = 2;
    ierr = KSPSetTolerances(sub_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,maxit);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"    ...ksp_ratio %g, new maxit %d\n\n",ksp_ratio,maxit);CHKERRQ(ierr);
  }
  check->its0 = its; /* save current outer KSP iteration number */
  PetscFunctionReturn(0);
}

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

  PetscFunctionBegin;
  if (spec->current) {
    ierr = KSPSolve(spec->kspcheap,ksp->vec_rhs,ksp->vec_sol);CHKERRQ(ierr);
    ierr = KSPSpecEstPropagateUp(ksp,spec->kspcheap);CHKERRQ(ierr);
  } else {
    PetscInt  i,its,neig;
    PetscReal *real,*imag,rad = 0;
    ierr = KSPSolve(spec->kspest,ksp->vec_rhs,ksp->vec_sol);CHKERRQ(ierr);
    ierr = KSPSpecEstPropagateUp(ksp,spec->kspest);CHKERRQ(ierr);
    ierr = KSPComputeExtremeSingularValues(spec->kspest,&spec->max,&spec->min);CHKERRQ(ierr);

    ierr = KSPGetIterationNumber(spec->kspest,&its);CHKERRQ(ierr);
    ierr = PetscMalloc2(its,PetscReal,&real,its,PetscReal,&imag);CHKERRQ(ierr);
    ierr = KSPComputeEigenvalues(spec->kspest,its,real,imag,&neig);CHKERRQ(ierr);
    for (i=0; i<neig; i++) {
      /* We would really like to compute w (nominally 1/radius) to minimize |1-wB|.  Empirically it
         is better to compute rad = |1-B| than rad = |B|.  There must be a cheap way to do better. */
      rad = PetscMax(rad,PetscRealPart(PetscSqrtScalar((PetscScalar)(PetscSqr(real[i]-1.) + PetscSqr(imag[i])))));
    }
    ierr = PetscFree2(real,imag);CHKERRQ(ierr);
    spec->radius = rad;

    ierr = KSPChebyshevSetEigenvalues(spec->kspcheap,spec->max*spec->maxfactor,spec->min*spec->minfactor);CHKERRQ(ierr);
    ierr = KSPRichardsonSetScale(spec->kspcheap,spec->richfactor/spec->radius);
    ierr = PetscInfo3(ksp,"Estimated singular value min=%G max=%G, spectral radius=%G",spec->min,spec->max,spec->radius);CHKERRQ(ierr);
    spec->current = PETSC_TRUE;
  }
  PetscFunctionReturn(0);
}

/**
   My handrolled SNES monitor.  Gives some information about KSP convergence as well
   as looking pretty.

   @param snes Petsc nonlinear context
   @param its number of iterations
   @param norm norm of nonlinear residual
   @param dctx application context
*/
PetscErrorCode MonitorFunction(SNES snes, PetscInt its, double norm, void *dctx)
{
  PetscErrorCode ierr;
  PetscInt lits;
  PetscMPIInt rank;
  KSP ksp;
  KSPConvergedReason kspreason;
  PetscReal kspnorm;

  PetscFunctionBegin;
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);

  ierr = SNESGetKSP(snes, &ksp);CHKERRQ(ierr);
  ierr = KSPGetConvergedReason(ksp, &kspreason);CHKERRQ(ierr);
  ierr = KSPGetResidualNorm(ksp, &kspnorm);CHKERRQ(ierr);
  ierr = KSPGetIterationNumber(ksp, &lits);CHKERRQ(ierr);

  ierr = PetscPrintf(PETSC_COMM_WORLD, "  %d SNES norm %e, %d KSP its last norm %e",
		       its, norm, lits, kspnorm);CHKERRQ(ierr);
  if (kspreason < 0) {
      ierr = PetscPrintf(PETSC_COMM_WORLD, ", KSP failed: %s", KSPConvergedReasons[kspreason]);CHKERRQ(ierr);
  }

  ierr = PetscPrintf(PETSC_COMM_WORLD, ".\n");CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

bool PotentialSolve::Solve(Vec delta, Vec pot, double bias) {
    PetscInt its;
    KSPConvergedReason reason;
    bool retval;

    // Delta to -Delta
    VecScale(delta, -1.0/bias);

    // Actually solve
    KSPSolve(solver, delta, pot);
    KSPGetConvergedReason(solver,&reason);
    if (reason<0) {
        PetscPrintf(PETSC_COMM_WORLD,"Diverged : %d.\n",reason);
        retval=false;
    } else {
        KSPGetIterationNumber(solver,&its);
        PetscPrintf(PETSC_COMM_WORLD,"\nConvergence in %d iterations.\n",(int)its);
        retval=true;
    }


    // Remove any values that might have crept in here
    _dg1.ZeroPad(pot);

    // Clean up
    VecScale(delta, -1.0*bias);

    return retval;
}

MGSolver_Status _GetSolveStatus( MGSolver_PETScData* mgData ) {
	PC			pc;
	const KSPType		kspType;
	const PCType		pcType;
	KSPConvergedReason	reason;
	PetscErrorCode		ec;

	ec = KSPGetType( mgData->ksp, &kspType );
	CheckPETScError( ec );
	ec = KSPGetPC( mgData->ksp, &pc );
	CheckPETScError( ec );
	ec = PCGetType( pc, &pcType );
	CheckPETScError( ec );

	if( !strcmp( kspType, KSPRICHARDSON ) && !strcmp( pcType, PCSOR ) ) {
		double		rnorm;
		PetscInt	curIt;

		//rnorm = PETScMatrixSolver_GetResidualNorm( self );
		//curIt = PETScMatrixSolver_GetIterations( self );
		rnorm = _GetResidualNorm( mgData );
		KSPGetIterationNumber( mgData->ksp, &curIt );
		//PETScMatrixSolver_SetNormType( self, MultigridSolver_NormType_Preconditioned );
		KSPSetNormType( mgData->ksp, MultigridSolver_NormType_Preconditioned );
		ec = KSPDefaultConverged( mgData->ksp, curIt, (PetscScalar)rnorm, &reason, PETSC_NULL );
		CheckPETScError( ec );
	}
	else {
		ec = KSPGetConvergedReason( mgData->ksp, &reason );
		CheckPETScError( ec );
	}

	return reason;
}

  /// Solve again w/ the specified RHS, put result in specified vector (specialized)
  void p_resolveImpl(const VectorType& b, VectorType& x) const
  {
    PetscErrorCode ierr(0);
    int me(this->processor_rank());
    try {
      const Vec *bvec(PETScVector(b));
      Vec *xvec(PETScVector(x));

      ierr = KSPSolve(p_KSP, *bvec, *xvec); CHKERRXX(ierr);
      int its;
      KSPConvergedReason reason;
      PetscReal rnorm;
      ierr = KSPGetIterationNumber(p_KSP, &its); CHKERRXX(ierr);
      ierr = KSPGetConvergedReason(p_KSP, &reason); CHKERRXX(ierr);
      ierr = KSPGetResidualNorm(p_KSP, &rnorm); CHKERRXX(ierr);
      std::string msg;
      if (reason < 0) {
        msg = 
          boost::str(boost::format("%d: PETSc KSP diverged after %d iterations, reason: %d") % 
                     me % its % reason);
        throw Exception(msg);
      } else if (me == 0) {
        msg = 
          boost::str(boost::format("%d: PETSc KSP converged after %d iterations, reason: %d") % 
                     me % its % reason);
        std::cerr << msg << std::endl;
      }
    } catch (const PETSC_EXCEPTION_TYPE& e) {
      throw PETScException(ierr, e);
    } catch (const Exception& e) {
      throw e;
    }
  }    

static PetscErrorCode TaoSolve_SSILS(Tao tao)
{
  TAO_SSLS                     *ssls = (TAO_SSLS *)tao->data;
  PetscReal                    psi, ndpsi, normd, innerd, t=0;
  PetscReal                    delta, rho;
  PetscInt                     iter=0,kspits;
  TaoConvergedReason           reason;
  TaoLineSearchConvergedReason ls_reason;
  PetscErrorCode               ierr;

  PetscFunctionBegin;
  /* Assume that Setup has been called!
     Set the structure for the Jacobian and create a linear solver. */
  delta = ssls->delta;
  rho = ssls->rho;

  ierr = TaoComputeVariableBounds(tao);CHKERRQ(ierr);
  ierr = VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);CHKERRQ(ierr);
  ierr = TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,Tao_SSLS_FunctionGradient,tao);CHKERRQ(ierr);
  ierr = TaoLineSearchSetObjectiveRoutine(tao->linesearch,Tao_SSLS_Function,tao);CHKERRQ(ierr);

  /* Calculate the function value and fischer function value at the
     current iterate */
  ierr = TaoLineSearchComputeObjectiveAndGradient(tao->linesearch,tao->solution,&psi,ssls->dpsi);CHKERRQ(ierr);
  ierr = VecNorm(ssls->dpsi,NORM_2,&ndpsi);CHKERRQ(ierr);

  while (1) {
    ierr=PetscInfo3(tao, "iter: %D, merit: %g, ndpsi: %g\n",iter, (double)ssls->merit, (double)ndpsi);CHKERRQ(ierr);
    /* Check the termination criteria */
    ierr = TaoMonitor(tao,iter++,ssls->merit,ndpsi,0.0,t,&reason);CHKERRQ(ierr);
    if (reason!=TAO_CONTINUE_ITERATING) break;

    /* Calculate direction.  (Really negative of newton direction.  Therefore,
       rest of the code uses -d.) */
    ierr = KSPSetOperators(tao->ksp,tao->jacobian,tao->jacobian_pre);CHKERRQ(ierr);
    ierr = KSPSolve(tao->ksp,ssls->ff,tao->stepdirection);CHKERRQ(ierr);
    ierr = KSPGetIterationNumber(tao->ksp,&kspits);CHKERRQ(ierr);
    tao->ksp_its+=kspits;
    ierr = VecNorm(tao->stepdirection,NORM_2,&normd);CHKERRQ(ierr);
    ierr = VecDot(tao->stepdirection,ssls->dpsi,&innerd);CHKERRQ(ierr);

    /* Make sure that we have a descent direction */
    if (innerd <= delta*pow(normd, rho)) {
      ierr = PetscInfo(tao, "newton direction not descent\n");CHKERRQ(ierr);
      ierr = VecCopy(ssls->dpsi,tao->stepdirection);CHKERRQ(ierr);
      ierr = VecDot(tao->stepdirection,ssls->dpsi,&innerd);CHKERRQ(ierr);
    }

    ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
    innerd = -innerd;

    ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
    ierr = TaoLineSearchApply(tao->linesearch,tao->solution,&psi,ssls->dpsi,tao->stepdirection,&t,&ls_reason);CHKERRQ(ierr);
    ierr = VecNorm(ssls->dpsi,NORM_2,&ndpsi);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

PetscErrorCode SNESQNApply_BadBroyden(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            *U  = qn->U;
  Vec            *T  = qn->V;

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

  PetscFunctionBegin;
  if (it < m) l = it;
  ierr = VecCopy(D,Y);CHKERRQ(ierr);
  if (l > 0) {
    k    = (it-1)%l;
    ierr = SNESLineSearchGetLambda(snes->linesearch,&qn->lambda[k]);CHKERRQ(ierr);
    ierr = VecCopy(Dold,U[k]);CHKERRQ(ierr);
    ierr = VecAXPY(U[k],-1.0,D);CHKERRQ(ierr);
    ierr = VecCopy(Xold,T[k]);CHKERRQ(ierr);
    ierr = VecAXPY(T[k],-1.0,X);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);
  }

  /* inward recursion starting at the first update and working forward */
  if (l > 0) {
    for (i = 0; i < l-1; i++) {
      j    = (it+i-l)%l;
      k    = (it+i-l+1)%l;
      h    = (it-1)%l;
      ierr = VecDotBegin(U[j],U[h],&gdot);CHKERRQ(ierr);
      ierr = VecDotBegin(U[j],U[j],&udot);CHKERRQ(ierr);
      ierr = VecDotEnd(U[j],U[h],&gdot);CHKERRQ(ierr);
      ierr = VecDotEnd(U[j],U[j],&udot);CHKERRQ(ierr);
      ierr = VecAXPY(Y,PetscRealPart(gdot)/PetscRealPart(udot),T[k]);CHKERRQ(ierr);
      ierr = VecAXPY(Y,-(1.-qn->lambda[k])*PetscRealPart(gdot)/PetscRealPart(udot),T[j]);CHKERRQ(ierr);
      if (qn->monitor) {
        ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
        ierr = PetscViewerASCIIPrintf(qn->monitor, "it: %d k: %d gdot: %14.12e\n", it, k, PetscRealPart(gdot));CHKERRQ(ierr);
        ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
      }
    }
  }
  PetscFunctionReturn(0);
}

static PetscErrorCode KSPSpecEstPropagateUp(KSP ksp,KSP subksp)
{
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = KSPGetConvergedReason(subksp,&ksp->reason);CHKERRQ(ierr);
  ierr = KSPGetIterationNumber(subksp,&ksp->its);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

static PetscErrorCode SNESSolve_KSPONLY(SNES snes)
{
  PetscErrorCode     ierr;
  PetscInt           lits;
  Vec                Y,X,F;

  PetscFunctionBegin;

  if (snes->xl || snes->xu || snes->ops->computevariablebounds) {
    SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);
  }

  snes->numFailures            = 0;
  snes->numLinearSolveFailures = 0;
  snes->reason                 = SNES_CONVERGED_ITERATING;
  snes->iter                   = 0;
  snes->norm                   = 0.0;

  X = snes->vec_sol;
  F = snes->vec_func;
  Y = snes->vec_sol_update;

  ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr);
  if (snes->numbermonitors) {
    PetscReal fnorm;
    ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr);
    ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr);
  }

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

  /* Solve J Y = F, where J is Jacobian matrix */
  ierr = SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr);
  ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr);
  ierr = KSPSolve(snes->ksp,F,Y);CHKERRQ(ierr);
  snes->reason = SNES_CONVERGED_ITS;
  SNESCheckKSPSolve(snes);

  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);
  snes->iter++;

  /* Take the computed step. */
  ierr = VecAXPY(X,-1.0,Y);CHKERRQ(ierr);
  if (snes->numbermonitors) {
    PetscReal fnorm;
    ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr);
    ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr);
    ierr = SNESMonitor(snes,1,fnorm);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

bool
PETScKrylovLinearSolver::solveSystem(SAMRAIVectorReal<NDIM, double>& x, SAMRAIVectorReal<NDIM, double>& b)
{
    IBTK_TIMER_START(t_solve_system);

#if !defined(NDEBUG)
    TBOX_ASSERT(d_A);
#endif
    int ierr;

    // Initialize the solver, when necessary.
    const bool deallocate_after_solve = !d_is_initialized;
    if (deallocate_after_solve) initializeSolverState(x, b);
#if !defined(NDEBUG)
    TBOX_ASSERT(d_petsc_ksp);
#endif
    resetKSPOptions();

    // Allocate scratch data.
    d_b->allocateVectorData();

    // Solve the system using a PETSc KSP object.
    d_b->copyVector(Pointer<SAMRAIVectorReal<NDIM, double> >(&b, false));
    d_A->setHomogeneousBc(d_homogeneous_bc);
    d_A->modifyRhsForBcs(*d_b);
    d_A->setHomogeneousBc(true);
    PETScSAMRAIVectorReal::replaceSAMRAIVector(d_petsc_x, Pointer<SAMRAIVectorReal<NDIM, double> >(&x, false));
    PETScSAMRAIVectorReal::replaceSAMRAIVector(d_petsc_b, d_b);
    ierr = KSPSolve(d_petsc_ksp, d_petsc_b, d_petsc_x);
    IBTK_CHKERRQ(ierr);
    d_A->setHomogeneousBc(d_homogeneous_bc);
    d_A->imposeSolBcs(x);

    // Get iterations count and residual norm.
    ierr = KSPGetIterationNumber(d_petsc_ksp, &d_current_iterations);
    IBTK_CHKERRQ(ierr);
    ierr = KSPGetResidualNorm(d_petsc_ksp, &d_current_residual_norm);
    IBTK_CHKERRQ(ierr);
    d_A->setHomogeneousBc(d_homogeneous_bc);

    // Determine the convergence reason.
    KSPConvergedReason reason;
    ierr = KSPGetConvergedReason(d_petsc_ksp, &reason);
    IBTK_CHKERRQ(ierr);
    const bool converged = (static_cast<int>(reason) > 0);
    if (d_enable_logging) reportKSPConvergedReason(reason, plog);

    // Dealocate scratch data.
    d_b->deallocateVectorData();

    // Deallocate the solver, when necessary.
    if (deallocate_after_solve) deallocateSolverState();

    IBTK_TIMER_STOP(t_solve_system);
    return converged;
} // solveSystem

PetscErrorCode SolveH(MPI_Comm comm, KSP ksp, Mat H, Vec rhs, Vec sol)
{
  /*-----------------KSP Solving------------------*/   
  PetscErrorCode ierr;
  PetscLogDouble t1,t2,tpast;
  ierr = Ptime(&t1);CHKERRQ(ierr);

  if (itsH>(maxit-5)){
    ierr = KSPSetOperators(ksp,H,H);CHKERRQ(ierr);}
  else{
    ierr = KSPSetReusePreconditioner(ksp,PETSC_TRUE);CHKERRQ(ierr);}

   ierr = KSPSolve(ksp,rhs,sol);CHKERRQ(ierr);
   ierr = KSPGetIterationNumber(ksp,&itsH);CHKERRQ(ierr);

   // if GMRES is stopped due to maxit, then redo it with sparse direct solve;
   if(itsH>(maxit-2))
     {
	ierr = KSPSetOperators(ksp,H,H);CHKERRQ(ierr);
	ierr = KSPSolve(ksp,rhs,sol);CHKERRQ(ierr);
	ierr = KSPGetIterationNumber(ksp,&itsH);CHKERRQ(ierr);
     }

  //Print kspsolving information
  double norm;
  Vec xdiff;
  ierr=VecDuplicate(sol,&xdiff);CHKERRQ(ierr);
  ierr = MatMult(H,sol, xdiff);CHKERRQ(ierr);
  ierr = VecAXPY(xdiff,-1.0,rhs);CHKERRQ(ierr);
  ierr = VecNorm(xdiff,NORM_INFINITY,&norm);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"==> Helmholtz filter solution: norm of error %g, Kryolv Iterations %d----\n ",norm,itsH);CHKERRQ(ierr);    

  ierr = Ptime(&t2);CHKERRQ(ierr);
  tpast = t2 - t1;

  int rank;
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  if(rank==0) PetscPrintf(PETSC_COMM_SELF,"==> Helmholtz filter solution: the runing time is %f s \n",tpast);
  /*--------------Finish KSP Solving---------------*/ 

  VecDestroy(&xdiff);
  PetscFunctionReturn(0);
}

PetscErrorCode BSSCR_PCBFBTSubKSPMonitor( KSP ksp, PetscInt index, PetscLogDouble time )
{
        PetscInt max_it;
        PetscReal rnorm;
        KSPConvergedReason reason;

        KSPGetIterationNumber( ksp, &max_it );
        KSPGetResidualNorm( ksp, &rnorm );
        KSPGetConvergedReason( ksp, &reason );
        PetscPrintf(((PetscObject)ksp)->comm,"    PCBFBTSubKSP (%d): %D Residual norm; r0 %12.12e, r %12.12e: Reason %s: Time %5.5e \n", 
		index, max_it, ksp->rnorm0, rnorm, KSPConvergedReasons[reason], time );

        PetscFunctionReturn(0);
}

int TaoLinearSolverPetsc::Solve(TaoVec* tv, TaoVec* tw, TaoTruth *flag)
{
  TaoVecPetsc *pv = dynamic_cast <TaoVecPetsc *> (tv);
  TaoVecPetsc *pw = dynamic_cast <TaoVecPetsc *> (tw);
  int info;
  PetscInt its;

  PetscFunctionBegin;

  info = KSPSolve(ksp,pv->GetVec(),pw->GetVec()); CHKERRQ(info);
  info = KSPGetIterationNumber(ksp, &its); CHKERRQ(info);

  this->linear_its=PetscMax(its,-its);
  if (its>0) *flag=TAO_TRUE; 
  else *flag=TAO_FALSE;

  PetscFunctionReturn(0);
}

static PetscErrorCode PCApply_KSP(PC pc,Vec x,Vec y)
{
  PetscErrorCode     ierr;
  PetscInt           its;
  PC_KSP             *jac = (PC_KSP*)pc->data;

  PetscFunctionBegin;
  if (jac->ksp->presolve) {
    ierr = VecCopy(x,y);CHKERRQ(ierr);
    ierr = KSPSolve(jac->ksp,y,y);CHKERRQ(ierr);
  } else {
    ierr = KSPSolve(jac->ksp,x,y);CHKERRQ(ierr);
  }
  ierr = KSPCheckSolve(jac->ksp,pc,y);CHKERRQ(ierr);
  ierr      = KSPGetIterationNumber(jac->ksp,&its);CHKERRQ(ierr);
  jac->its += its;
  PetscFunctionReturn(0);
}

PetscErrorCode BSSCR_Lp_monitor( KSP ksp, PetscInt index )
{
	PetscInt max_it;
	PetscReal rnorm;
	KSPConvergedReason reason;
	
	KSPGetIterationNumber( ksp, &max_it );
	KSPGetResidualNorm( ksp, &rnorm );
	KSPGetConvergedReason( ksp, &reason );
	if (ksp->reason > 0) {
		PetscPrintf(((PetscObject)ksp)->comm,"\t<Lp(%d)>: Linear solve converged. its.=%.4d ; |r|=%5.5e ; Reason=%s\n", 
				index, max_it, rnorm, KSPConvergedReasons[reason] );
	} else {
		PetscPrintf(((PetscObject)ksp)->comm,"\t<Lp(%d)>: Linear solve did not converge. its.=%.4d ; |r|=%5.5e ; Reason=%s\n", 
				index, max_it, rnorm, KSPConvergedReasons[reason]);
	}
	
	PetscFunctionReturn(0);
}

/*
 * Must be passed a KSP solver that has "converged", with KSPSetComputeEigenvalues() called before the solve
 */
static PetscErrorCode KSPChebyshevComputeExtremeEigenvalues_Private(KSP kspest,PetscReal *emin,PetscReal *emax)
{
  PetscErrorCode ierr;
  PetscInt       n,neig;
  PetscReal      *re,*im,min,max;

  PetscFunctionBegin;
  ierr = KSPGetIterationNumber(kspest,&n);CHKERRQ(ierr);
  ierr = PetscMalloc2(n,&re,n,&im);CHKERRQ(ierr);
  ierr = KSPComputeEigenvalues(kspest,n,re,im,&neig);CHKERRQ(ierr);
  min  = PETSC_MAX_REAL;
  max  = PETSC_MIN_REAL;
  for (n=0; n<neig; n++) {
    min = PetscMin(min,re[n]);
    max = PetscMax(max,re[n]);
  }
  ierr  = PetscFree2(re,im);CHKERRQ(ierr);
  *emax = max;
  *emin = min;
  PetscFunctionReturn(0);
}

//.........这里部分代码省略.........
  ierr = KSPSetOperators(cksp,cmat,cmat);CHKERRQ(ierr);
  ierr = KSPGetPC(cksp,&pc);CHKERRQ(ierr);
  ierr = PCSetType(pc,PCLU);CHKERRQ(ierr);
  ierr = KSPSetType(cksp,KSPPREONLY);CHKERRQ(ierr);

  /* zero is finest level */
  for (i=0; i<levels-1; i++) {
    ierr = PCMGSetResidual(pcmg,levels - 1 - i,residual,(Mat)0);CHKERRQ(ierr);
    ierr = MatCreateShell(PETSC_COMM_WORLD,N[i+1],N[i],N[i+1],N[i],(void*)0,&mat[i]);CHKERRQ(ierr);
    ierr = MatShellSetOperation(mat[i],MATOP_MULT,(void (*)(void))restrct);CHKERRQ(ierr);
    ierr = MatShellSetOperation(mat[i],MATOP_MULT_TRANSPOSE_ADD,(void (*)(void))interpolate);CHKERRQ(ierr);
    ierr = PCMGSetInterpolation(pcmg,levels - 1 - i,mat[i]);CHKERRQ(ierr);
    ierr = PCMGSetRestriction(pcmg,levels - 1 - i,mat[i]);CHKERRQ(ierr);
    ierr = PCMGSetCyclesOnLevel(pcmg,levels - 1 - i,cycles);CHKERRQ(ierr);

    /* set smoother */
    ierr = PCMGGetSmoother(pcmg,levels - 1 - i,&ksp[i]);CHKERRQ(ierr);
    ierr = KSPGetPC(ksp[i],&pc);CHKERRQ(ierr);
    ierr = PCSetType(pc,PCSHELL);CHKERRQ(ierr);
    ierr = PCShellSetName(pc,"user_precond");CHKERRQ(ierr);
    ierr = PCShellGetName(pc,&shellname);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"level=%D, PCShell name is %s\n",i,shellname);CHKERRQ(ierr);

    /* this is a dummy! since KSP requires a matrix passed in  */
    ierr = KSPSetOperators(ksp[i],mat[i],mat[i]);CHKERRQ(ierr);
    /*
        We override the matrix passed in by forcing it to use Richardson with
        a user provided application. This is non-standard and this practice
        should be avoided.
    */
    ierr = PCShellSetApplyRichardson(pc,gauss_seidel);CHKERRQ(ierr);
    if (use_jacobi) {
      ierr = PCShellSetApplyRichardson(pc,jacobi);CHKERRQ(ierr);
    }
    ierr = KSPSetType(ksp[i],KSPRICHARDSON);CHKERRQ(ierr);
    ierr = KSPSetInitialGuessNonzero(ksp[i],PETSC_TRUE);CHKERRQ(ierr);
    ierr = KSPSetTolerances(ksp[i],PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,smooths);CHKERRQ(ierr);

    ierr = VecCreateSeq(PETSC_COMM_SELF,N[i],&x);CHKERRQ(ierr);

    X[levels - 1 - i] = x;
    if (i > 0) {
      ierr = PCMGSetX(pcmg,levels - 1 - i,x);CHKERRQ(ierr);
    }
    ierr = VecCreateSeq(PETSC_COMM_SELF,N[i],&x);CHKERRQ(ierr);

    B[levels -1 - i] = x;
    if (i > 0) {
      ierr = PCMGSetRhs(pcmg,levels - 1 - i,x);CHKERRQ(ierr);
    }
    ierr = VecCreateSeq(PETSC_COMM_SELF,N[i],&x);CHKERRQ(ierr);

    R[levels - 1 - i] = x;

    ierr = PCMGSetR(pcmg,levels - 1 - i,x);CHKERRQ(ierr);
  }
  /* create coarse level vectors */
  ierr = VecCreateSeq(PETSC_COMM_SELF,N[levels-1],&x);CHKERRQ(ierr);
  ierr = PCMGSetX(pcmg,0,x);CHKERRQ(ierr); X[0] = x;
  ierr = VecCreateSeq(PETSC_COMM_SELF,N[levels-1],&x);CHKERRQ(ierr);
  ierr = PCMGSetRhs(pcmg,0,x);CHKERRQ(ierr); B[0] = x;

  /* create matrix multiply for finest level */
  ierr = MatCreateShell(PETSC_COMM_WORLD,N[0],N[0],N[0],N[0],(void*)0,&fmat);CHKERRQ(ierr);
  ierr = MatShellSetOperation(fmat,MATOP_MULT,(void (*)(void))amult);CHKERRQ(ierr);
  ierr = KSPSetOperators(kspmg,fmat,fmat);CHKERRQ(ierr);

  ierr = CalculateSolution(N[0],&solution);CHKERRQ(ierr);
  ierr = CalculateRhs(B[levels-1]);CHKERRQ(ierr);
  ierr = VecSet(X[levels-1],0.0);CHKERRQ(ierr);

  ierr = residual((Mat)0,B[levels-1],X[levels-1],R[levels-1]);CHKERRQ(ierr);
  ierr = CalculateError(solution,X[levels-1],R[levels-1],e);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_SELF,"l_2 error %g max error %g resi %g\n",(double)e[0],(double)e[1],(double)e[2]);CHKERRQ(ierr);

  ierr = KSPSolve(kspmg,B[levels-1],X[levels-1]);CHKERRQ(ierr);
  ierr = KSPGetIterationNumber(kspmg,&its);CHKERRQ(ierr);
  ierr = residual((Mat)0,B[levels-1],X[levels-1],R[levels-1]);CHKERRQ(ierr);
  ierr = CalculateError(solution,X[levels-1],R[levels-1],e);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_SELF,"its %D l_2 error %g max error %g resi %g\n",its,(double)e[0],(double)e[1],(double)e[2]);CHKERRQ(ierr);

  ierr = PetscFree(N);CHKERRQ(ierr);
  ierr = VecDestroy(&solution);CHKERRQ(ierr);

  /* note we have to keep a list of all vectors allocated, this is
     not ideal, but putting it in MGDestroy is not so good either*/
  for (i=0; i<levels; i++) {
    ierr = VecDestroy(&X[i]);CHKERRQ(ierr);
    ierr = VecDestroy(&B[i]);CHKERRQ(ierr);
    if (i) {ierr = VecDestroy(&R[i]);CHKERRQ(ierr);}
  }
  for (i=0; i<levels-1; i++) {
    ierr = MatDestroy(&mat[i]);CHKERRQ(ierr);
  }
  ierr = MatDestroy(&cmat);CHKERRQ(ierr);
  ierr = MatDestroy(&fmat);CHKERRQ(ierr);
  ierr = KSPDestroy(&kspmg);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return 0;
}

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

  /*
     Create parallel vectors compatible with the DMDA.
  */
  ierr = DMCreateGlobalVector(da,&u);CHKERRQ(ierr);
  ierr = VecDuplicate(u,&b);CHKERRQ(ierr);
  ierr = VecDuplicate(u,&x);CHKERRQ(ierr);

  /*
     Set exact solution; then compute right-hand-side vector.
     By default we use an exact solution of a vector with all
     elements of 1.0;  Alternatively, using the runtime option
     -random_sol forms a solution vector with random components.
  */
  ierr = PetscOptionsGetBool(NULL,"-random_exact_sol",&flg,NULL);CHKERRQ(ierr);
  if (flg) {
    ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);CHKERRQ(ierr);
    ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
    ierr = VecSetRandom(u,rctx);CHKERRQ(ierr);
    ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr);
  } else {
    ierr = VecSet(u,1.);CHKERRQ(ierr);
  }
  ierr = MatMult(A,u,b);CHKERRQ(ierr);

  /*
     View the exact solution vector if desired
  */
  flg  = PETSC_FALSE;
  ierr = PetscOptionsGetBool(NULL,"-view_exact_sol",&flg,NULL);CHKERRQ(ierr);
  if (flg) {ierr = VecView(u,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);}

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                Create the linear solver and set various options
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  /*
     Create linear solver context
  */
  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);

  /*
     Set operators. Here the matrix that defines the linear system
     also serves as the preconditioning matrix.
  */
  ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);

  /*
    Set runtime options, e.g.,
        -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
    These options will override those specified above as long as
    KSPSetFromOptions() is called _after_ any other customization
    routines.
  */
  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                      Solve the linear system
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                      Check solution and clean up
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  /*
     Check the error
  */
  ierr = VecAXPY(x,-1.,u);CHKERRQ(ierr);
  ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr);
  ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);

  /*
     Print convergence information.  PetscPrintf() produces a single
     print statement from all processes that share a communicator.
     An alternative is PetscFPrintf(), which prints to a file.
  */
  ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm,its);CHKERRQ(ierr);

  /*
     Free work space.  All PETSc objects should be destroyed when they
     are no longer needed.
  */
  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
  ierr = VecDestroy(&u);CHKERRQ(ierr);
  ierr = VecDestroy(&x);CHKERRQ(ierr);
  ierr = VecDestroy(&b);CHKERRQ(ierr);
  ierr = MatDestroy(&A);CHKERRQ(ierr);
  ierr = DMDestroy(&da);CHKERRQ(ierr);

  /*
     Always call PetscFinalize() before exiting a program.  This routine
       - finalizes the PETSc libraries as well as MPI
       - provides summary and diagnostic information if certain runtime
         options are chosen (e.g., -log_summary).
  */
  ierr = PetscFinalize();
  return 0;
}

int main(int argc,char **args)
{
  Vec            x,b,u;      /* approx solution, RHS, exact solution */
  Mat            A;            /* linear system matrix */
  KSP            ksp;         /* KSP context */
  PetscErrorCode ierr;
  PetscInt       n = 10,its, dim,p = 1,use_random;
  PetscScalar    none = -1.0,pfive = 0.5;
  PetscReal      norm;
  PetscRandom    rctx;
  TestType       type;
  PetscBool      flg;

  PetscInitialize(&argc,&args,(char *)0,help);
  ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetInt(PETSC_NULL,"-p",&p,PETSC_NULL);CHKERRQ(ierr);
  switch (p) {
    case 1:  type = TEST_1;      dim = n;   break;
    case 2:  type = TEST_2;      dim = n;   break;
    case 3:  type = TEST_3;      dim = n;   break;
    case 4:  type = HELMHOLTZ_1; dim = n*n; break;
    case 5:  type = HELMHOLTZ_2; dim = n*n; break;
    default: type = TEST_1;      dim = n;
  }

  /* Create vectors */
  ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
  ierr = VecSetSizes(x,PETSC_DECIDE,dim);CHKERRQ(ierr);
  ierr = VecSetFromOptions(x);CHKERRQ(ierr);
  ierr = VecDuplicate(x,&b);CHKERRQ(ierr);
  ierr = VecDuplicate(x,&u);CHKERRQ(ierr);

  use_random = 1;
  flg        = PETSC_FALSE;
  ierr = PetscOptionsGetBool(PETSC_NULL,"-norandom",&flg,PETSC_NULL);CHKERRQ(ierr);
  if (flg) {
    use_random = 0;
    ierr = VecSet(u,pfive);CHKERRQ(ierr);
  } else {
    ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);CHKERRQ(ierr);
    ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
    ierr = VecSetRandom(u,rctx);CHKERRQ(ierr);
  }

  /* Create and assemble matrix */
  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,dim,dim);CHKERRQ(ierr);
  ierr = MatSetFromOptions(A);CHKERRQ(ierr);
  ierr = MatSetUp(A);CHKERRQ(ierr);
  ierr = FormTestMatrix(A,n,type);CHKERRQ(ierr);
  ierr = MatMult(A,u,b);CHKERRQ(ierr);
  flg  = PETSC_FALSE;
  ierr = PetscOptionsGetBool(PETSC_NULL,"-printout",&flg,PETSC_NULL);CHKERRQ(ierr);
  if (flg) {
    ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
    ierr = VecView(u,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
    ierr = VecView(b,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  }

  /* Create KSP context; set operators and options; solve linear system */
  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
  ierr = KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
  ierr = KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

  /* Check error */
  ierr = VecAXPY(x,none,u);CHKERRQ(ierr);
  ierr  = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr);
  ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %G,Iterations %D\n",norm,its);CHKERRQ(ierr);

  /* Free work space */
  ierr = VecDestroy(&x);CHKERRQ(ierr); ierr = VecDestroy(&u);CHKERRQ(ierr);
  ierr = VecDestroy(&b);CHKERRQ(ierr); ierr = MatDestroy(&A);CHKERRQ(ierr);
  if (use_random) {ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr);}
  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return 0;
}