C++ Ifpack类代码示例

// ===================================================
// Methods
// ===================================================
Int
PreconditionerIfpack::buildPreconditioner( operator_type& matrix )
{
    M_operator     = matrix->matrixPtr();

    M_overlapLevel = this->M_list.get( "overlap level", -1 );

    M_precType     = this->M_list.get( "prectype", "Amesos" );

    Ifpack factory;

    M_preconditioner.reset( factory.Create( M_precType, M_operator.get(), M_overlapLevel ) );

    M_precType += "_Ifpack";

    if ( !M_preconditioner.get() )
    {
        ERROR_MSG( "Preconditioner not set, something went wrong in its computation\n" );
    }

    IFPACK_CHK_ERR( M_preconditioner->SetParameters( this->M_list ) );
    IFPACK_CHK_ERR( M_preconditioner->Initialize() );
    IFPACK_CHK_ERR( M_preconditioner->Compute() );

    this->M_preconditionerCreated = true;

    return ( EXIT_SUCCESS );
}

  void IfpackSmoother::Setup(Level &currentLevel) {
    FactoryMonitor m(*this, "Setup Smoother", currentLevel);
    if (SmootherPrototype::IsSetup() == true)
      GetOStream(Warnings0, 0) << "Warning: MueLu::IfpackSmoother::Setup(): Setup() has already been called";

    A_ = Factory::Get< RCP<Matrix> >(currentLevel, "A");

    double lambdaMax = -1.0;
    if (type_ == "Chebyshev") {
      std::string maxEigString   = "chebyshev: max eigenvalue";
      std::string eigRatioString = "chebyshev: ratio eigenvalue";

      try {
        lambdaMax = Teuchos::getValue<Scalar>(this->GetParameter(maxEigString));
        this->GetOStream(Statistics1, 0) << maxEigString << " (cached with smoother parameter list) = " << lambdaMax << std::endl;

      } catch (Teuchos::Exceptions::InvalidParameterName) {
        lambdaMax = A_->GetMaxEigenvalueEstimate();

        if (lambdaMax != -1.0) {
          this->GetOStream(Statistics1, 0) << maxEigString << " (cached with matrix) = " << lambdaMax << std::endl;
          this->SetParameter(maxEigString, ParameterEntry(lambdaMax));
        }
      }

      // Calculate the eigenvalue ratio
      const Scalar defaultEigRatio = 20;

      Scalar ratio = defaultEigRatio;
      try {
        ratio = Teuchos::getValue<Scalar>(this->GetParameter(eigRatioString));

      } catch (Teuchos::Exceptions::InvalidParameterName) {
        this->SetParameter(eigRatioString, ParameterEntry(ratio));
      }

      if (currentLevel.GetLevelID()) {
        // Update ratio to be
        //   ratio = max(number of fine DOFs / number of coarse DOFs, defaultValue)
        //
        // NOTE: We don't need to request previous level matrix as we know for sure it was constructed
        RCP<const Matrix> fineA = currentLevel.GetPreviousLevel()->Get<RCP<Matrix> >("A");
        size_t nRowsFine   = fineA->getGlobalNumRows();
        size_t nRowsCoarse = A_->getGlobalNumRows();

        ratio = std::max(ratio, as<Scalar>(nRowsFine)/nRowsCoarse);

        this->GetOStream(Statistics1, 0) << eigRatioString << " (computed) = " << ratio << std::endl;
        this->SetParameter(eigRatioString, ParameterEntry(ratio));
      }
    }

    RCP<Epetra_CrsMatrix> epA = Utils::Op2NonConstEpetraCrs(A_);

    Ifpack factory;
    prec_ = rcp(factory.Create(type_, &(*epA), overlap_));
    TEUCHOS_TEST_FOR_EXCEPTION(prec_.is_null(), Exceptions::RuntimeError, "Could not create an Ifpack preconditioner with type = \"" << type_ << "\"");
    SetPrecParameters();
    prec_->Compute();

    SmootherPrototype::IsSetup(true);

    if (type_ == "Chebyshev" && lambdaMax == -1.0) {
      Teuchos::RCP<Ifpack_Chebyshev> chebyPrec = rcp_dynamic_cast<Ifpack_Chebyshev>(prec_);
      if (chebyPrec != Teuchos::null) {
        lambdaMax = chebyPrec->GetLambdaMax();
        A_->SetMaxEigenvalueEstimate(lambdaMax);
        this->GetOStream(Statistics1, 0) << "chebyshev: max eigenvalue (calculated by Ifpack)" << " = " << lambdaMax << std::endl;
      }
      TEUCHOS_TEST_FOR_EXCEPTION(lambdaMax == -1.0, Exceptions::RuntimeError, "MueLu::IfpackSmoother::Setup(): no maximum eigenvalue estimate");
    }

    this->GetOStream(Statistics0, 0) << description() << std::endl;
  }

int main(int argc, char *argv[]) {
  //
  int MyPID = 0;
#ifdef EPETRA_MPI
  // Initialize MPI
  MPI_Init(&argc,&argv);
  Epetra_MpiComm Comm(MPI_COMM_WORLD);
  MyPID = Comm.MyPID();
#else
  Epetra_SerialComm Comm;
#endif
  //
  typedef double                            ST;
  typedef Teuchos::ScalarTraits<ST>        SCT;
  typedef SCT::magnitudeType                MT;
  typedef Epetra_MultiVector                MV;
  typedef Epetra_Operator                   OP;
  typedef Belos::MultiVecTraits<ST,MV>     MVT;
  typedef Belos::OperatorTraits<ST,MV,OP>  OPT;

  using Teuchos::ParameterList;
  using Teuchos::RCP;
  using Teuchos::rcp;

  bool success = false;
  bool verbose = false;
  try {
    bool proc_verbose = false;
    int frequency = -1;        // frequency of status test output.
    int blocksize = 1;         // blocksize
    int numrhs = 1;            // number of right-hand sides to solve for
    int maxrestarts = 15;      // maximum number of restarts allowed
    int maxiters = -1;         // maximum number of iterations allowed per linear system
    int maxsubspace = 25;      // maximum number of blocks the solver can use for the subspace
    std::string filename("orsirr1.hb");
    MT tol = 1.0e-5;           // relative residual tolerance

    // Specify whether to use RHS as initial guess. If false, use zero.
    bool useRHSAsInitialGuess = false;

    Teuchos::CommandLineProcessor cmdp(false,true);
    cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
    cmdp.setOption("use-rhs","use-zero",&useRHSAsInitialGuess,"Use RHS as initial guess.");
    cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters).");
    cmdp.setOption("filename",&filename,"Filename for test matrix.  Acceptable file extensions: *.hb,*.mtx,*.triU,*.triS");
    cmdp.setOption("tol",&tol,"Relative residual tolerance used by GMRES solver.");
    cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
    cmdp.setOption("block-size",&blocksize,"Block size used by GMRES.");
    cmdp.setOption("max-iters",&maxiters,"Maximum number of iterations per linear system (-1 = adapted to problem/block size).");
    cmdp.setOption("max-subspace",&maxsubspace,"Maximum number of blocks the solver can use for the subspace.");
    cmdp.setOption("max-restarts",&maxrestarts,"Maximum number of restarts allowed for GMRES solver.");
    if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
      return EXIT_FAILURE;
    }

    if (!verbose)
      frequency = -1;  // reset frequency if test is not verbose

    //
    // *************Get the problem*********************
    //
    RCP<Epetra_Map> Map;
    RCP<Epetra_CrsMatrix> A;
    RCP<Epetra_MultiVector> B, X;
    RCP<Epetra_Vector> vecB, vecX;
    EpetraExt::readEpetraLinearSystem(filename, Comm, &A, &Map, &vecX, &vecB);
    A->OptimizeStorage();
    proc_verbose = verbose && (MyPID==0);  /* Only print on the zero processor */

    // Check to see if the number of right-hand sides is the same as requested.
    if (numrhs>1) {
      X = rcp( new Epetra_MultiVector( *Map, numrhs ) );
      B = rcp( new Epetra_MultiVector( *Map, numrhs ) );
      X->Seed();
      X->Random();
      OPT::Apply( *A, *X, *B );
      X->PutScalar( 0.0 );
    }
    else {
      X = Teuchos::rcp_implicit_cast<Epetra_MultiVector>(vecX);
      B = Teuchos::rcp_implicit_cast<Epetra_MultiVector>(vecB);
    }

    // If requested, use a copy of B as initial guess
    if (useRHSAsInitialGuess)
    {
      X->Update(1.0, *B, 0.0);
    }

    //
    // ************Construct preconditioner*************
    //
    ParameterList ifpackList;

    // allocates an IFPACK factory. No data is associated
    // to this object (only method Create()).
    Ifpack Factory;

    // create the preconditioner. For valid PrecType values,
    // please check the documentation
//.........这里部分代码省略.........

int main(int argc, char *argv[]) {
    //
    int MyPID = 0;
#ifdef EPETRA_MPI
    // Initialize MPI
    MPI_Init(&argc,&argv);
    Epetra_MpiComm Comm(MPI_COMM_WORLD);
    MyPID = Comm.MyPID();
#else
    Epetra_SerialComm Comm;
#endif
    //
    typedef double                            ST;
    typedef Teuchos::ScalarTraits<ST>        SCT;
    typedef SCT::magnitudeType                MT;
    typedef Epetra_MultiVector                MV;
    typedef Epetra_Operator                   OP;
    typedef Belos::MultiVecTraits<ST,MV>     MVT;
    typedef Belos::OperatorTraits<ST,MV,OP>  OPT;

    using Teuchos::ParameterList;
    using Teuchos::RCP;
    using Teuchos::rcp;

    bool verbose = false, debug = false, proc_verbose = false, strict_conv = false;
    int frequency = -1;        // frequency of status test output.
    int blocksize = 1;         // blocksize
    int numrhs = 1;            // number of right-hand sides to solve for
    int maxiters = -1;         // maximum number of iterations allowed per linear system
    int maxsubspace = 50;      // maximum number of blocks the solver can use for the subspace
    int maxrestarts = 15;      // number of restarts allowed
    std::string filename("orsirr1.hb");
    std::string precond("none");
    MT tol = 1.0e-5;           // relative residual tolerance
    MT polytol = tol/10;       // relative residual tolerance for polynomial construction

    Teuchos::CommandLineProcessor cmdp(false,true);
    cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
    cmdp.setOption("debug","nondebug",&debug,"Print debugging information from solver.");
    cmdp.setOption("strict-conv","not-strict-conv",&strict_conv,"Require solver to strictly adhere to convergence tolerance.");
    cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters).");
    cmdp.setOption("filename",&filename,"Filename for test matrix.  Acceptable file extensions: *.hb,*.mtx,*.triU,*.triS");
    cmdp.setOption("precond",&precond,"Preconditioning type (none, left, right).");
    cmdp.setOption("tol",&tol,"Relative residual tolerance used by GMRES solver.");
    cmdp.setOption("poly-tol",&polytol,"Relative residual tolerance used to construct the GMRES polynomial.");
    cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
    cmdp.setOption("block-size",&blocksize,"Block size used by GMRES.");
    cmdp.setOption("max-iters",&maxiters,"Maximum number of iterations per linear system (-1 = adapted to problem/block size).");
    cmdp.setOption("max-subspace",&maxsubspace,"Maximum number of blocks the solver can use for the subspace.");
    cmdp.setOption("max-restarts",&maxrestarts,"Maximum number of restarts allowed for GMRES solver.");
    if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
        return -1;
    }
    if (!verbose)
        frequency = -1;  // reset frequency if test is not verbose
    //
    // Get the problem
    //
    RCP<Epetra_Map> Map;
    RCP<Epetra_CrsMatrix> A;
    RCP<Epetra_MultiVector> B, X;
    RCP<Epetra_Vector> vecB, vecX;
    EpetraExt::readEpetraLinearSystem(filename, Comm, &A, &Map, &vecX, &vecB);
    A->OptimizeStorage();
    proc_verbose = verbose && (MyPID==0);  /* Only print on the zero processor */

    // Check to see if the number of right-hand sides is the same as requested.
    if (numrhs>1) {
        X = rcp( new Epetra_MultiVector( *Map, numrhs ) );
        B = rcp( new Epetra_MultiVector( *Map, numrhs ) );
        X->Seed();
        X->Random();
        OPT::Apply( *A, *X, *B );
        X->PutScalar( 0.0 );
    }
    else {
        X = Teuchos::rcp_implicit_cast<Epetra_MultiVector>(vecX);
        B = Teuchos::rcp_implicit_cast<Epetra_MultiVector>(vecB);
    }
    //
    // ************Construct preconditioner*************
    //
    RCP<Belos::EpetraPrecOp> belosPrec;

    if (precond != "none") {
        ParameterList ifpackList;

        // allocates an IFPACK factory. No data is associated
        // to this object (only method Create()).
        Ifpack Factory;

        // create the preconditioner. For valid PrecType values,
        // please check the documentation
        std::string PrecType = "ILU"; // incomplete LU
        int OverlapLevel = 1; // must be >= 0. If Comm.NumProc() == 1,
        // it is ignored.

        RCP<Ifpack_Preconditioner> Prec = Teuchos::rcp( Factory.Create(PrecType, &*A, OverlapLevel) );
        assert(Prec != Teuchos::null);

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

//.........这里部分代码省略.........
      double value = (double) ( globNumCol -1 - gid );
      int numEntries = 1;
      vecX->ReplaceGlobalValues( numEntries, &value, &gid );
    }
    bool Trans = false;
    A->Multiply( Trans, *vecX, *vecB ); // Create a consistent linear system
    // At this point, the initial guess is exact.
    bool zeroInitGuess = false; // annihilate initial guess
    bool goodInitGuess = true; // initial guess near solution
    if( zeroInitGuess )
      {
        vecX->PutScalar( 0.0 );
      }
    else
      {
        if( goodInitGuess )
          {
            double value = 1.e-2; // "Rel RHS Err" and "Rel Mat Err" apply to the residual equation,
            int numEntries = 1;   // norm( b - A x_k ) ?<? relResTol norm( b- Axo).
            int index = 0;        // norm(b) is inaccessible to LSQR.
            vecX->SumIntoMyValues(  numEntries, &value, &index);
          }
      }
    X = Teuchos::rcp_implicit_cast<Epetra_MultiVector>(vecX);
    B = Teuchos::rcp_implicit_cast<Epetra_MultiVector>(vecB);
  }
  //
  // ************Construct preconditioner*************
  //
  ParameterList ifpackList;

  // allocates an IFPACK factory. No data is associated
  // to this object (only method Create()).
  Ifpack Factory;  // do support transpose multiply

  // create the preconditioner. For valid PrecType values,
  // please check the documentation
  std::string PrecType = "ILU"; // incomplete LU
  int OverlapLevel = 1; // nonnegative

  RCP<Ifpack_Preconditioner> Prec = Teuchos::rcp( Factory.Create(PrecType, &*A, OverlapLevel) );
  assert(Prec != Teuchos::null);

  // specify parameters for ILU
  ifpackList.set("fact: level-of-fill", 1);
  // the combine mode is on the following:
  // "Add", "Zero", "Insert", "InsertAdd", "Average", "AbsMax"
  // Their meaning is as defined in file Epetra_CombineMode.h
  ifpackList.set("schwarz: combine mode", "Add");
  // sets the parameters
  IFPACK_CHK_ERR(Prec->SetParameters(ifpackList));

  // initialize the preconditioner. At this point the matrix must
  // have been FillComplete()'d, but actual values are ignored.
  IFPACK_CHK_ERR(Prec->Initialize());

  // Builds the preconditioners, by looking for the values of
  // the matrix.
  IFPACK_CHK_ERR(Prec->Compute());

  {
    const int errcode = Prec->SetUseTranspose (true);
    if (errcode != 0) {
      throw std::logic_error ("Oh hai! Ifpack_Preconditioner doesn't know how to apply its transpose.");
    } else {
      (void) Prec->SetUseTranspose (false);

int main(int argc, char *argv[])
{

  // Initialize MPI
#ifdef HAVE_MPI
  MPI_Init(&argc,&argv);
#endif

  // Create a communicator for Epetra objects
#ifdef HAVE_MPI
  Epetra_MpiComm Comm( MPI_COMM_WORLD );
#else
  Epetra_SerialComm Comm;
#endif

  // Get the process ID and the total number of processors
  int MyPID = Comm.MyPID();
  int NumProc = Comm.NumProc();

  // Check verbosity level
  bool verbose = false;
  if (argc > 1)
    if (argv[1][0]=='-' && argv[1][1]=='v')
      verbose = true;

  // Get the number of elements from the command line
  int NumGlobalElements = 0;
  if ((argc > 2) && (verbose))
    NumGlobalElements = atoi(argv[2]) + 1;
  else if ((argc > 1) && (!verbose))
    NumGlobalElements = atoi(argv[1]) + 1;
  else
    NumGlobalElements = 101;

  // The number of unknowns must be at least equal to the
  // number of processors.
  if (NumGlobalElements < NumProc) {
    std::cout << "numGlobalBlocks = " << NumGlobalElements
      << " cannot be < number of processors = " << NumProc << std::endl;
    std::cout << "Test failed!" << std::endl;
    throw "NOX Error";
  }

  bool success = false;

  try {
    // Create the interface between NOX and the application
    // This object is derived from NOX::Epetra::Interface
    Teuchos::RCP<Interface> interface =
      Teuchos::rcp(new Interface(NumGlobalElements, Comm));

    // Set the PDE factor (for nonlinear forcing term).  This could be specified
    // via user input.
    interface->setPDEfactor(1000.0);

    // Begin Nonlinear Solver ************************************

    // Create the top level parameter list
    Teuchos::RCP<Teuchos::ParameterList> IfpackParamsPtr =
      Teuchos::rcp(new Teuchos::ParameterList);

    // Set the printing parameters in the "Printing" sublist
    Teuchos::ParameterList printParams;
    printParams.set("MyPID", MyPID);
    printParams.set("Output Precision", 3);
    printParams.set("Output Processor", 0);
    if (verbose)
      printParams.set("Output Information",
          NOX::Utils::OuterIteration +
          NOX::Utils::OuterIterationStatusTest +
          NOX::Utils::InnerIteration +
          NOX::Utils::LinearSolverDetails +
          NOX::Utils::Parameters +
          NOX::Utils::Details +
          NOX::Utils::Warning +
          NOX::Utils::Debug +
          NOX::Utils::TestDetails +
          NOX::Utils::Error);
    else
      printParams.set("Output Information", NOX::Utils::Error +
          NOX::Utils::TestDetails);

    // Create a print class for controlling output below
    NOX::Utils p(printParams);


    // *******************************
    // Setup Test Objects
    // *******************************

    // Create Linear Objects
    // Get the vector from the Problem
    if (verbose)
      p.out() << "Creating Vectors and Matrices" << std::endl;
    Teuchos::RCP<Epetra_Vector> solution_vec =
      interface->getSolution();
    Teuchos::RCP<Epetra_Vector> rhs_vec =
      Teuchos::rcp(new Epetra_Vector(*solution_vec));
    Teuchos::RCP<Epetra_Vector> lhs_vec =
      Teuchos::rcp(new Epetra_Vector(*solution_vec));
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
#ifdef HAVE_MPI
  MPI_Init(&argc,&argv);
  Epetra_MpiComm comm(MPI_COMM_WORLD);
#else
  Epetra_SerialComm comm;
#endif

  Galeri::core::Workspace::setNumDimensions(3);

  Galeri::grid::Loadable domain, boundary;

  int numGlobalElementsX = 2 * comm.NumProc();
  int numGlobalElementsY = 2;
  int numGlobalElementsZ = 2;

  int mx = comm.NumProc();
  int my = 1;
  int mz = 1;

  Galeri::grid::Generator::
  getCubeWithHexs(comm, numGlobalElementsX, numGlobalElementsY, numGlobalElementsZ,
                  mx, my, mz, domain, boundary);

  Epetra_Map matrixMap(domain.getNumGlobalVertices(), 0, comm);

  Epetra_FECrsMatrix A(Copy, matrixMap, 0);
  Epetra_FEVector    LHS(matrixMap);
  Epetra_FEVector    RHS(matrixMap);

  Galeri::problem::ScalarLaplacian<Laplacian> problem("Hex", 1, 8);

  problem.integrate(domain, A, RHS);

  LHS.PutScalar(0.0);

  problem.imposeDirichletBoundaryConditions(boundary, A, RHS, LHS);

  // ============================================================ //
  // Solving the linear system is the next step, using the IFPACK //
  // factory. This is done by using the IFPACK factory, then      //
  // asking for IC preconditioner, and setting few parameters     //
  // using a Teuchos::ParameterList.                              //
  // ============================================================ //
  
  Ifpack Factory;
  Ifpack_Preconditioner* Prec = Factory.Create("IC", &A, 0);

  Teuchos::ParameterList list;
  
  list.set("fact: level-of-fill", 1);
  IFPACK_CHK_ERR(Prec->SetParameters(list));
  IFPACK_CHK_ERR(Prec->Initialize());
  IFPACK_CHK_ERR(Prec->Compute());

  Epetra_LinearProblem linearProblem(&A, &LHS, &RHS);

  AztecOO solver(linearProblem);
  solver.SetAztecOption(AZ_solver, AZ_cg);
  solver.SetPrecOperator(Prec);
  solver.Iterate(1550, 1e-9);

  // visualization using MEDIT -- a VTK module is available as well
  Galeri::viz::MEDIT::write(domain, "sol", LHS);

  // now compute the norm of the solution
  problem.computeNorms(domain, LHS);

#ifdef HAVE_MPI
  MPI_Finalize();
#endif
}

int Ifpack_IHSS::Compute(){
  if(!IsInitialized_) Initialize();

  int rv;
  Ifpack Factory;
  Epetra_CrsMatrix *Askew=0,*Aherm=0;
  Ifpack_Preconditioner *Pskew=0, *Pherm=0;
  Time_.ResetStartTime();

  // Create Aherm (w/o diagonal)
  rv=EpetraExt::MatrixMatrix::Add(*A_,false,.5,*A_,true,.5,Aherm);
  Aherm->FillComplete();
  if(rv) IFPACK_CHK_ERR(-1); 

  // Grab Aherm's diagonal
  Epetra_Vector avec(Aherm->RowMap());
  IFPACK_CHK_ERR(Aherm->ExtractDiagonalCopy(avec));


  // Compute alpha using the Bai, Golub & Ng 2003 formula, not the more multigrid-appropriate Hamilton, Benzi and Haber 2007.
  //  PowerMethod(Aherm, EigMaxIters_,LambdaMax_);
  //  Alpha_=LambdaMax_ / sqrt(EigRatio_);

  // Try something more Hamilton inspired, using the maximum diagonal value of Aherm.
  avec.MaxValue(&Alpha_);

  // Add alpha to the diagonal of Aherm
  for(int i=0;i<Aherm->NumMyRows();i++) avec[i]+=Alpha_;   
  IFPACK_CHK_ERR(Aherm->ReplaceDiagonalValues(avec));
  Aherm_=rcp(Aherm);

  // Compute Askew (and add diagonal)
  Askew=new Epetra_CrsMatrix(Copy,A_->RowMap(),0);
  rv=EpetraExt::MatrixMatrix::Add(*A_,false,.5,*A_,true,-.5,Askew);
  if(rv) IFPACK_CHK_ERR(-2); 
  for(int i=0;i<Askew->NumMyRows();i++) {
    int gid=Askew->GRID(i);
    Askew->InsertGlobalValues(gid,1,&Alpha_,&gid);
  }
  Askew->FillComplete();
  Askew_=rcp(Askew);

  // Compute preconditioner for Aherm
  Teuchos::ParameterList PLh=List_.sublist("ihss: hermetian list");
  string htype=List_.get("ihss: hermetian type","ILU");
  Pherm= Factory.Create(htype, Aherm);
  Pherm->SetParameters(PLh);
  IFPACK_CHK_ERR(Pherm->Compute());
  Pherm_=rcp(Pherm);

  // Compute preconditoner for Askew 
  Teuchos::ParameterList PLs=List_.sublist("ihss: skew hermetian list");
  string stype=List_.get("ihss: skew hermetian type","ILU");
  Pskew= Factory.Create(stype, Askew);
  Pskew->SetParameters(PLs);
  IFPACK_CHK_ERR(Pskew->Compute());
  Pskew_=rcp(Pskew);
  
  // Label
  sprintf(Label_, "IFPACK IHSS (H,S)=(%s/%s)",htype.c_str(),stype.c_str()); 

  // Counters
  IsComputed_=true;
  NumCompute_++;
  ComputeTime_ += Time_.ElapsedTime();
  return 0;
}

int main(int argc, char *argv[]) {
  //
#ifdef EPETRA_MPI	
  // Initialize MPI	
  MPI_Init(&argc,&argv); 	
  Belos::MPIFinalize mpiFinalize; // Will call finalize with *any* return
  (void)mpiFinalize;
#endif
  //
  using Teuchos::RCP;
  using Teuchos::rcp;
  //
  // Get test parameters from command-line processor
  //  
  bool verbose = false, proc_verbose = false;
  int frequency = -1; // how often residuals are printed by solver 
  int numrhs = 15;  // total number of right-hand sides to solve for
  int blocksize = 10;  // blocksize used by solver
  int maxiters = -1; // maximum number of iterations for the solver to use
  std::string filename("bcsstk14.hb");
  double tol = 1.0e-5;  // relative residual tolerance

  Teuchos::CommandLineProcessor cmdp(false,true);
  cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
  cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters).");
  cmdp.setOption("filename",&filename,"Filename for Harwell-Boeing test matrix.");
  cmdp.setOption("tol",&tol,"Relative residual tolerance used by CG solver.");
  cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
  cmdp.setOption("block-size",&blocksize,"Block size to be used by CG solver.");
  cmdp.setOption("max-iters",&maxiters,"Maximum number of iterations per linear system (-1 := adapted to problem/block size).");

  if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
    return -1;
  }
  if (!verbose)
    frequency = -1;  // Reset frequency if verbosity is off
  //
  // Get the problem
  //
  int MyPID;
  RCP<Epetra_CrsMatrix> A;
  RCP<Epetra_MultiVector> X, B;
  int return_val =Belos::createEpetraProblem(filename,NULL,&A,&B,&X,&MyPID);
  if(return_val != 0) return return_val;
  proc_verbose = ( verbose && (MyPID==0) );
  //
  // Solve using Belos
  //
  typedef double                           ST;
  typedef Epetra_Operator                  OP;
  typedef Epetra_MultiVector               MV;
  typedef Belos::OperatorTraits<ST,MV,OP> OPT;
  typedef Belos::MultiVecTraits<ST,MV>    MVT;
  //
  // *****Construct initial guess and random right-hand-sides *****
  //
  if (numrhs != 1) {
    X = rcp( new Epetra_MultiVector( A->Map(), numrhs ) );
    MVT::MvRandom( *X );
    B = rcp( new Epetra_MultiVector( A->Map(), numrhs ) );
    OPT::Apply( *A, *X, *B );
    MVT::MvInit( *X, 0.0 );
  }
  //
  // ************Construct preconditioner*************
  //
  ParameterList ifpackList;

  // allocates an IFPACK factory. No data is associated
  // to this object (only method Create()).
  Ifpack Factory;

  // create the preconditioner. For valid PrecType values,
  // please check the documentation
  std::string PrecType = "ICT"; // incomplete Cholesky 
  int OverlapLevel = 0; // must be >= 0. If Comm.NumProc() == 1,
                        // it is ignored.

  RCP<Ifpack_Preconditioner> Prec = Teuchos::rcp( Factory.Create(PrecType, &*A, OverlapLevel) );
  assert(Prec != Teuchos::null);

  // specify parameters for ICT
  ifpackList.set("fact: drop tolerance", 1e-9);
  ifpackList.set("fact: ict level-of-fill", 1.0);
  // the combine mode is on the following:
  // "Add", "Zero", "Insert", "InsertAdd", "Average", "AbsMax"
  // Their meaning is as defined in file Epetra_CombineMode.h
  ifpackList.set("schwarz: combine mode", "Add");
  // sets the parameters
  IFPACK_CHK_ERR(Prec->SetParameters(ifpackList));

  // initialize the preconditioner. At this point the matrix must
  // have been FillComplete()'d, but actual values are ignored.
  IFPACK_CHK_ERR(Prec->Initialize());

  // Builds the preconditioners, by looking for the values of
  // the matrix.
  IFPACK_CHK_ERR(Prec->Compute());

  // Create the Belos preconditioned operator from the Ifpack preconditioner.
//.........这里部分代码省略.........

void MxGeoMultigridPrec::setup() {

  std::string smoothType = pList.get("linear solver : smoother type", "Gauss-Seidel");
  std::string coarseSmoothType = pList.get("linear solver : coarse smoother", "Amesos");
  int smoothSweeps = pList.get("linear solver : smoother sweeps", 1);
  int overlap = 1;
  
  levels = pList.get("linear solver : levels", 1);
  cycles = pList.get("linear solver : cycles", 1);
  output = pList.get("linear solver : output", 0);
  rcf = pList.get("linear solver : remove const field", false);

  // Gauss-Seidel preconditioner options
  Teuchos::ParameterList GSList;
  GSList.set("relaxation: type", "Gauss-Seidel");
  //GSList.set("relaxation: type", "Jacobi");
  GSList.set("relaxation: sweeps", smoothSweeps);
  GSList.set("relaxation: damping factor", 1.);
  GSList.set("relaxation: zero starting solution", false);

  // Chebyshev preconditioner options
  Teuchos::ParameterList ChebList;
  ChebList.set("chebyshev: degree", smoothSweeps);
  ChebList.set("chebyshev: zero starting solution", false);
  ChebList.set("chebyshev: ratio eigenvalue", 30.);

  // Amesos direct solver parameters
  Teuchos::ParameterList AmesosList;
  AmesosList.set("amesos: solver type", "Amesos_Klu");
  AmesosList.set("AddToDiag", 1.e-12);



  // smoother parameters for intermediate levels
  Teuchos::ParameterList vSmootherList;
  std::string vPrecType;
  if (smoothType == "Gauss-Seidel") {
    vPrecType = "point relaxation";
    vSmootherList = GSList;
  }
  else if (smoothType == "Chebyshev") {
    vPrecType = smoothType;
    vSmootherList = ChebList;
  }
  else {
    std::cout << "MxGeoMultigridPrec::setup(): unknown smoother type, '" << smoothType << "' for intermediate smoothers. Using gauss-seidel.\n";
    vSmootherList = GSList;
  }

  // smoother parameters for coarsest level
  Teuchos::ParameterList cSmootherList;
  std::string cPrecType;
  if (coarseSmoothType == "Gauss-Seidel") {
    cSmootherList = GSList;
    cPrecType = "point relaxation";
  }
  else if (coarseSmoothType == "Chebyshev") {
    cSmootherList = ChebList;
    cPrecType = coarseSmoothType;
  }
  else if (coarseSmoothType == "Amesos") {
    cSmootherList = AmesosList;
    cPrecType = coarseSmoothType;
  }
  else {
    std::cout << "MxGeoMultigridPrec::setup(): unknown smoother type, '" << coarseSmoothType << "', for coarsest smoother. Using Amesos_Klu.\n";
    cSmootherList = AmesosList;
    cPrecType = "Amesos";
  }

  Teuchos::ParameterList smootherList;
  std::string precType;

  Ifpack factory;
  
  double maxEig, minEig;

  // special case for when Amesos wants to alter our original operator
  if (levels == 1 and coarseSmoothType == "Amesos") {
    fineOpCopy = Teuchos::rcp(new Epetra_CrsMatrix(*ops[0])); //apparently Amesos changes the input matrix?
    smoothers.push_back(Teuchos::rcp(new Ifpack_Amesos(&*fineOpCopy)));
    smoothers[0]->SetParameters(AmesosList);
    smoothers[0]->Initialize();
    smoothers[0]->Compute();
  }
  else {
  //level 0 is original operator
    for (int level = 0; level < levels; ++level) {
      std::cout << "  Initializing multigrid level " << level << std::endl;

      if (level == levels - 1) {
        smootherList = cSmootherList;
        precType = cPrecType;
      }
      else {
        smootherList = vSmootherList;
        precType = vPrecType;
      }

      // always find max eigenvalue for at least the prolongator smoother,
//.........这里部分代码省略.........

int main(int argc, char *argv[]) {

#ifdef HAVE_MPI
  MPI_Init(&argc,&argv);
  Epetra_MpiComm Comm (MPI_COMM_WORLD);
#else
  Epetra_SerialComm Comm;
#endif

  Teuchos::ParameterList GaleriList;
  int nx = 30; 

  GaleriList.set("nx", nx);
  //  GaleriList.set("ny", nx * Comm.NumProc());
  GaleriList.set("ny", nx);
  GaleriList.set("mx", 1);
  GaleriList.set("my", Comm.NumProc());
  GaleriList.set("alpha", .0);
  GaleriList.set("diff", 1.0);
  GaleriList.set("conv", 100.0);

  Teuchos::RefCountPtr<Epetra_Map> Map = Teuchos::rcp( Galeri::CreateMap64("Cartesian2D", Comm, GaleriList) );
  Teuchos::RefCountPtr<Epetra_CrsMatrix> A = Teuchos::rcp( Galeri::CreateCrsMatrix("UniFlow2D", &*Map, GaleriList) );
  Teuchos::RefCountPtr<Epetra_MultiVector> LHS = Teuchos::rcp( new Epetra_MultiVector(*Map, 1) );
  Teuchos::RefCountPtr<Epetra_MultiVector> RHS = Teuchos::rcp( new Epetra_MultiVector(*Map, 1) );
  LHS->PutScalar(0.0); RHS->Random();
  Ifpack Factory;  
  int Niters = 100;

  // ============================= //
  // Construct IHSS preconditioner //
  // ============================= //
  Teuchos::RefCountPtr<Ifpack_Preconditioner> Prec = Teuchos::rcp( Factory.Create("IHSS", &*A,0) );
  Teuchos::ParameterList List;
  List.set("ihss: hermetian type","ILU");
  List.set("ihss: skew hermetian type","ILU");
  List.set("ihss: ratio eigenvalue",100.0);
  // Could set sublist values here to better control the ILU, but this isn't needed for this example.
  IFPACK_CHK_ERR(Prec->SetParameters(List));
  IFPACK_CHK_ERR(Prec->Compute());

  // ============================= //
  // Create solver Object          //
  // ============================= //

  AztecOO solver;
  solver.SetUserMatrix(&*A);
  solver.SetLHS(&*LHS);
  solver.SetRHS(&*RHS);
  solver.SetAztecOption(AZ_solver,AZ_gmres);
  solver.SetPrecOperator(&*Prec);
  solver.SetAztecOption(AZ_output, 1); 
  solver.Iterate(Niters, 1e-8);

  // ============================= //
  // Construct SORa preconditioner //
  // ============================= //
  Teuchos::RefCountPtr<Ifpack_Preconditioner> Prec2 = Teuchos::rcp( Factory.Create("SORa", &*A,0) );
  Teuchos::ParameterList List2;
  List2.set("sora: sweeps",1);
  // Could set sublist values here to better control the ILU, but this isn't needed for this example.
  IFPACK_CHK_ERR(Prec2->SetParameters(List2));
  IFPACK_CHK_ERR(Prec2->Compute());

  // ============================= //
  // Create solver Object          //
  // ============================= //
  AztecOO solver2;
  LHS->PutScalar(0.0);
  solver2.SetUserMatrix(&*A);
  solver2.SetLHS(&*LHS);
  solver2.SetRHS(&*RHS);
  solver2.SetAztecOption(AZ_solver,AZ_gmres);
  solver2.SetPrecOperator(&*Prec2);
  solver2.SetAztecOption(AZ_output, 1); 
  solver2.Iterate(Niters, 1e-8);

#ifdef HAVE_MPI
  MPI_Finalize() ;
#endif

  return(EXIT_SUCCESS);
}

int main(int argc, char *argv[]) {
  using std::cout;
  using std::endl;
  //
  bool haveM = true;

#ifdef EPETRA_MPI
  // Initialize MPI
  MPI_Init(&argc,&argv);
  Epetra_MpiComm Comm( MPI_COMM_WORLD );
#else
  Epetra_SerialComm Comm;
#endif

  int MyPID = Comm.MyPID();

  bool verbose=false;
  bool isHermitian=false;
  std::string k_filename = "bfw782a.mtx";
  std::string m_filename = "bfw782b.mtx";
  std::string which = "LR";
  Teuchos::CommandLineProcessor cmdp(false,true);
  cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
  cmdp.setOption("sort",&which,"Targetted eigenvalues (SM,LM,SR,or LR).");
  cmdp.setOption("K-filename",&k_filename,"Filename and path of the stiffness matrix.");
  cmdp.setOption("M-filename",&m_filename,"Filename and path of the mass matrix.");
  if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
#ifdef HAVE_MPI
    MPI_Finalize();
#endif
    return -1;
  }
  //
  //**********************************************************************
  //******************Set up the problem to be solved*********************
  //**********************************************************************
  //
  // *****Read in matrix from file******
  //
  Teuchos::RCP<Epetra_Map> Map;
  Teuchos::RCP<Epetra_CrsMatrix> K, M;
  EpetraExt::readEpetraLinearSystem( k_filename, Comm, &K, &Map );

  if (haveM) {
    EpetraExt::readEpetraLinearSystem( m_filename, Comm, &M, &Map );
  }
  //
  // Build Preconditioner
  //
  Ifpack factory;
  std::string ifpack_type = "ILUT";
  int overlap = 0;
  Teuchos::RCP<Ifpack_Preconditioner> ifpack_prec = Teuchos::rcp(
          factory.Create( ifpack_type, K.get(), overlap ) );

  //
  // Set parameters and compute preconditioner
  //
  Teuchos::ParameterList ifpack_params;
  double droptol = 1e-2;
  double fill = 2.0;
  ifpack_params.set("fact: drop tolerance",droptol);
  ifpack_params.set("fact: ilut level-of-fill",fill);
  ifpack_prec->SetParameters(ifpack_params);
  ifpack_prec->Initialize();
  ifpack_prec->Compute();

  //
  // GeneralizedDavidson expects preconditioner to be applied with
  //  "Apply" rather than "Apply_Inverse"
  //
  Teuchos::RCP<Epetra_Operator> prec = Teuchos::rcp(
          new Epetra_InvOperator(ifpack_prec.get()) );

  //
  //************************************
  // Start the block Davidson iteration
  //***********************************
  //
  //  Variables used for the Block Arnoldi Method
  //
  int nev = 5;
  int blockSize = 5;
  int maxDim = 40;
  int maxRestarts = 10;
  double tol = 1.0e-8;

  // Set verbosity level
  int verbosity = Anasazi::Errors + Anasazi::Warnings;
  if (verbose) {
    verbosity += Anasazi::FinalSummary + Anasazi::TimingDetails;
  }
  //
  // Create parameter list to pass into solver
  //
  Teuchos::ParameterList MyPL;
  MyPL.set( "Verbosity", verbosity );
  MyPL.set( "Which", which );
  MyPL.set( "Block Size", blockSize );
  MyPL.set( "Maximum Subspace Dimension", maxDim );
//.........这里部分代码省略.........

int main(int argc, char *argv[]) {

#ifdef HAVE_MPI
  MPI_Init(&argc,&argv);
  Epetra_MpiComm Comm (MPI_COMM_WORLD);
#else
  Epetra_SerialComm Comm;
#endif

  int MyPID = Comm.MyPID();
  bool verbose = false; 
  if (MyPID==0) verbose = true;

  // The problem is defined on a 2D grid, global size is nx * nx.
  int nx = 30;
  Teuchos::ParameterList GaleriList;
  GaleriList.set("nx", nx);
  GaleriList.set("ny", nx * Comm.NumProc());
  GaleriList.set("mx", 1);
  GaleriList.set("my", Comm.NumProc());
  Teuchos::RefCountPtr<Epetra_Map> Map = Teuchos::rcp( Galeri::CreateMap("Cartesian2D", Comm, GaleriList) );
  Teuchos::RefCountPtr<Epetra_CrsMatrix> A = Teuchos::rcp( Galeri::CreateCrsMatrix("Laplace2D", &*Map, GaleriList) );
  Teuchos::RefCountPtr<Epetra_MultiVector> LHS = Teuchos::rcp( new Epetra_MultiVector(*Map, 1) );
  Teuchos::RefCountPtr<Epetra_MultiVector> RHS = Teuchos::rcp( new Epetra_MultiVector(*Map, 1) );
  LHS->PutScalar(0.0); RHS->Random();

  // ============================ //
  // Construct ILU preconditioner //
  // ---------------------------- //

  // I wanna test funky values to be sure that they have the same
  // influence on the algorithms, both old and new
  int    LevelFill = 2;
  double DropTol = 0.3333;
  double Condest;
  
  Teuchos::RefCountPtr<Ifpack_CrsIct> ICT;
  ICT = Teuchos::rcp( new Ifpack_CrsIct(*A,DropTol,LevelFill) );
  ICT->SetAbsoluteThreshold(0.00123);
  ICT->SetRelativeThreshold(0.9876);
  // Init values from A
  ICT->InitValues(*A);
  // compute the factors
  ICT->Factor();
  // and now estimate the condition number
  ICT->Condest(false,Condest);
  
  if( Comm.MyPID() == 0 ) {
    cout << "Condition number estimate (level-of-fill = "
	 << LevelFill <<  ") = " << Condest << endl;
  }

  // Define label for printing out during the solve phase
  string label = "Ifpack_CrsIct Preconditioner: LevelFill = " + toString(LevelFill) + 
                                                 " Overlap = 0"; 
  ICT->SetLabel(label.c_str());
  
  // Here we create an AztecOO object
  LHS->PutScalar(0.0);

  int Niters = 1200;

  AztecOO solver;
  solver.SetUserMatrix(&*A);
  solver.SetLHS(&*LHS);
  solver.SetRHS(&*RHS);
  solver.SetAztecOption(AZ_solver,AZ_cg);
  solver.SetPrecOperator(&*ICT);
  solver.SetAztecOption(AZ_output, 16); 
  solver.Iterate(Niters, 5.0e-5);

  int OldIters = solver.NumIters();

  // now rebuild the same preconditioner using ICT, we expect the same
  // number of iterations

  Ifpack Factory;
  Teuchos::RefCountPtr<Ifpack_Preconditioner> Prec = Teuchos::rcp( Factory.Create("IC", &*A) );

  Teuchos::ParameterList List;
  List.get("fact: level-of-fill", 2);
  List.get("fact: drop tolerance", 0.3333);
  List.get("fact: absolute threshold", 0.00123);
  List.get("fact: relative threshold", 0.9876);
  List.get("fact: relaxation value", 0.0);

  IFPACK_CHK_ERR(Prec->SetParameters(List));
  IFPACK_CHK_ERR(Prec->Compute());

  // Here we create an AztecOO object
  LHS->PutScalar(0.0);

  solver.SetUserMatrix(&*A);
  solver.SetLHS(&*LHS);
  solver.SetRHS(&*RHS);
  solver.SetAztecOption(AZ_solver,AZ_cg);
  solver.SetPrecOperator(&*Prec);
  solver.SetAztecOption(AZ_output, 16); 
  solver.Iterate(Niters, 5.0e-5);

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

//.........这里部分代码省略.........
                        Jac_vbr->BeginReplaceMyValues(global_row, 1, &global_col);
                        if (global_row==global_col)
                            Jac_vbr->SubmitBlockEntry(diag_block_matrix);
                        else
                            Jac_vbr->SubmitBlockEntry(block_matrix);
                        Jac_vbr->EndSubmitEntries();
                    }
                }
            }
            Jac_vbr->FillComplete();
            x = rcp(new Epetra_Vector(map));
            delta_x = rcp(new Epetra_Vector(map));
            f = rcp(new Epetra_Vector(map));

            x->PutScalar(1.0);
            Jac_vbr->Apply(*x,*f);

            Jac =
                rcpWithEmbeddedObjPostDestroy(new Epetra_VbrRowMatrix(Jac_vbr.get()),
                                              Jac_vbr);

        }

        if (print_debug_info) {
            x->Print(cout);
            Jac->Print(cout);
            f->Print(cout);
        }

        // *********************************************************
        // * Build Preconditioner (Ifpack)
        // *********************************************************

        Ifpack Factory;
        std::string PrecType = "ILU";
        int OverlapLevel = 1;
        RCP<Ifpack_Preconditioner> Prec =
            Teuchos::rcp( Factory.Create(PrecType, &*Jac, OverlapLevel) );
        ParameterList ifpackList;
        ifpackList.set("fact: drop tolerance", 1e-9);
        ifpackList.set("fact: level-of-fill", 1);
        ifpackList.set("schwarz: combine mode", "Add");
        IFPACK_CHK_ERR(Prec->SetParameters(ifpackList));
        IFPACK_CHK_ERR(Prec->Initialize());
        IFPACK_CHK_ERR(Prec->Compute());
        IFPACK_CHK_ERR(Prec->Condest());
        RCP<Belos::EpetraPrecOp> belosPrec =
            rcp( new Belos::EpetraPrecOp( Prec ) );

        // *********************************************************
        // * Build linear solver (Belos)
        // *********************************************************

        // Linear solver parameters
        typedef double                            ST;
        typedef Teuchos::ScalarTraits<ST>        SCT;
        typedef SCT::magnitudeType                MT;
        typedef Epetra_MultiVector                MV;
        typedef Epetra_Operator                   OP;
        typedef Belos::MultiVecTraits<ST,MV>     MVT;
        typedef Belos::OperatorTraits<ST,MV,OP>  OPT;

        RCP<ParameterList> belosList = rcp(new ParameterList);
        belosList->set<int>("Num Blocks", num_dof);
        belosList->set<int>("Block Size", 1);
        belosList->set<int>("Maximum Iterations", 400);

int main(int argc, char *argv[])
{

  // initialize MPI and Epetra communicator
#ifdef HAVE_MPI
  MPI_Init(&argc,&argv);
  Epetra_MpiComm Comm( MPI_COMM_WORLD );
#else
  Epetra_SerialComm Comm;
#endif

  Teuchos::ParameterList GaleriList;

  // The problem is defined on a 2D grid, global size is nx * nx.
  int nx = 30; 
  GaleriList.set("nx", nx);
  GaleriList.set("ny", nx * Comm.NumProc());
  GaleriList.set("mx", 1);
  GaleriList.set("my", Comm.NumProc());
  Teuchos::RefCountPtr<Epetra_Map> Map = Teuchos::rcp( Galeri::CreateMap("Cartesian2D", Comm, GaleriList) );
  Teuchos::RefCountPtr<Epetra_RowMatrix> A = Teuchos::rcp( Galeri::CreateCrsMatrix("Laplace2D", &*Map, GaleriList) );

  // =============================================================== //
  // B E G I N N I N G   O F   I F P A C K   C O N S T R U C T I O N //
  // =============================================================== //

  Teuchos::ParameterList List;

  // allocates an IFPACK factory. No data is associated 
  // to this object (only method Create()).
  Ifpack Factory;

  // create the preconditioner. For valid PrecType values,
  // please check the documentation
  std::string PrecType = "Amesos";
  int OverlapLevel = 2; // must be >= 0. If Comm.NumProc() == 1,
                        // it is ignored.

  Teuchos::RefCountPtr<Ifpack_Preconditioner> Prec = Teuchos::rcp( Factory.Create(PrecType, &*A, OverlapLevel) );
  assert(Prec != Teuchos::null);

  // specify the Amesos solver to be used. 
  // If the selected solver is not available,
  // IFPACK will try to use Amesos' KLU (which is usually always
  // compiled). Amesos' serial solvers are:
  // "Amesos_Klu", "Amesos_Umfpack", "Amesos_Superlu"
  List.set("amesos: solver type", "Amesos_Klu");

  // sets the parameters
  IFPACK_CHK_ERR(Prec->SetParameters(List));

  // initialize the preconditioner. At this point the matrix must
  // have been FillComplete()'d, but actual values are ignored.
  // At this call, Amesos will perform the symbolic factorization.
  IFPACK_CHK_ERR(Prec->Initialize());

  // Builds the preconditioners, by looking for the values of 
  // the matrix. At this call, Amesos will perform the
  // numeric factorization.
  IFPACK_CHK_ERR(Prec->Compute());

  // =================================================== //
  // E N D   O F   I F P A C K   C O N S T R U C T I O N //
  // =================================================== //

  // At this point, we need some additional objects
  // to define and solve the linear system.

  // defines LHS and RHS
  Epetra_Vector LHS(A->OperatorDomainMap());
  Epetra_Vector RHS(A->OperatorDomainMap());

  // solution is constant
  LHS.PutScalar(1.0);
  // now build corresponding RHS
  A->Apply(LHS,RHS);

  // now randomize the solution
  RHS.Random();

  // need an Epetra_LinearProblem to define AztecOO solver
  Epetra_LinearProblem Problem(&*A,&LHS,&RHS);

  // now we can allocate the AztecOO solver
  AztecOO Solver(Problem);

  // specify solver
  Solver.SetAztecOption(AZ_solver,AZ_gmres);
  Solver.SetAztecOption(AZ_output,32);

  // HERE WE SET THE IFPACK PRECONDITIONER
  Solver.SetPrecOperator(&*Prec);

  // .. and here we solve
  // NOTE: with one process, the solver must converge in
  // one iteration.
  Solver.Iterate(1550,1e-8);

#ifdef HAVE_MPI
  MPI_Finalize() ; 
//.........这里部分代码省略.........