C++ Epetra_LinearProblem::SetOperator方法代码示例

// ====================================================================== 
void Eig(const Operator& Op, MultiVector& ER, MultiVector& EI)
{
  int ierr;
  if (Op.GetDomainSpace() != Op.GetRangeSpace())
    ML_THROW("Matrix is not square", -1);

  ER.Reshape(Op.GetDomainSpace());
  EI.Reshape(Op.GetDomainSpace());

  Epetra_LinearProblem Problem;
  Problem.SetOperator(const_cast<Epetra_RowMatrix*>(Op.GetRowMatrix()));
  Amesos_Lapack Lapack(Problem);

  Epetra_Vector ER_Epetra(Op.GetRowMatrix()->RowMatrixRowMap());
  Epetra_Vector EI_Epetra(Op.GetRowMatrix()->RowMatrixRowMap());

  ierr = Lapack.GEEV(ER_Epetra, EI_Epetra);

  if (ierr)
    ML_THROW("GEEV returned error code = " + GetString(ierr), -1);
  
  for (int i = 0 ; i < ER.GetMyLength() ; ++i) {
    ER(i) = ER_Epetra[i];
    EI(i) = EI_Epetra[i];
  }
}

// ======================================================================
void Krylov(const Operator& A, const MultiVector& LHS,
            const MultiVector& RHS, const BaseOperator& Prec,
            Teuchos::ParameterList& List)
{
#ifndef HAVE_ML_AZTECOO
      std::cerr << "Please configure ML with --enable-aztecoo to use" << std::endl;
      std::cerr << "MLAPI Krylov solvers" << std::endl;
      exit(EXIT_FAILURE);
#else
  if (LHS.GetNumVectors() != 1)
    ML_THROW("FIXME: only one vector is currently supported", -1);

  Epetra_LinearProblem Problem;

  const Epetra_RowMatrix& A_Epetra = *(A.GetRowMatrix());

  Epetra_Vector LHS_Epetra(View,A_Epetra.OperatorDomainMap(),
                           (double*)&(LHS(0)));
  Epetra_Vector RHS_Epetra(View,A_Epetra.OperatorRangeMap(),
                           (double*)&(RHS(0)));

  // FIXME: this works only for Epetra-based operators
  Problem.SetOperator((const_cast<Epetra_RowMatrix*>(&A_Epetra)));
  Problem.SetLHS(&LHS_Epetra);
  Problem.SetRHS(&RHS_Epetra);

  AztecOO solver(Problem);

  EpetraBaseOperator Prec_Epetra(A_Epetra.OperatorDomainMap(),Prec);
  solver.SetPrecOperator(&Prec_Epetra);

  // get options from List
  int    NumIters = List.get("krylov: max iterations", 1550);
  double Tol      = List.get("krylov: tolerance", 1e-9);
  std::string type     = List.get("krylov: type", "gmres");
  int    output   = List.get("krylov: output level", GetPrintLevel());

  // set options in `solver'
  if (type == "cg")
    solver.SetAztecOption(AZ_solver, AZ_cg);
  else if (type == "cg_condnum")
    solver.SetAztecOption(AZ_solver, AZ_cg_condnum);
  else if (type == "gmres")
    solver.SetAztecOption(AZ_solver, AZ_gmres);
  else if (type == "gmres_condnum")
    solver.SetAztecOption(AZ_solver, AZ_gmres_condnum);
  else if (type == "fixed point")
    solver.SetAztecOption(AZ_solver, AZ_fixed_pt);
  else
    ML_THROW("krylov: type has incorrect value (" +
             type + ")", -1);

  solver.SetAztecOption(AZ_output, output);
  solver.Iterate(NumIters, Tol);
#endif

}

int PartialFactorizationOneStep( const char* AmesosClass,
				 const Epetra_Comm &Comm, 
				 bool transpose, 
				 bool verbose, 
				 Teuchos::ParameterList ParamList, 
				 Epetra_CrsMatrix *& Amat, 
				 double Rcond, 
				 int Steps ) 
{
	
  assert( Steps >= 0 && Steps < MaxNumSteps ) ; 

  int iam = Comm.MyPID() ; 
  int errors = 0 ; 

  const Epetra_Map *RangeMap = 
    transpose?&Amat->OperatorDomainMap():&Amat->OperatorRangeMap() ; 
  const Epetra_Map *DomainMap = 
    transpose?&Amat->OperatorRangeMap():&Amat->OperatorDomainMap() ; 

  Epetra_Vector xexact(*DomainMap);
  Epetra_Vector x(*DomainMap);

  Epetra_Vector b(*RangeMap);
  Epetra_Vector bcheck(*RangeMap);

  Epetra_Vector difference(*DomainMap);

  Epetra_LinearProblem Problem;
  Amesos_BaseSolver* Abase ; 
  Amesos Afactory;

  Abase = Afactory.Create( AmesosClass, Problem ) ; 

  std::string AC = AmesosClass ;
  if ( AC == "Amesos_Mumps" ) { 
    ParamList.set( "NoDestroy", true );
   Abase->SetParameters( ParamList ) ; 
  }

  double relresidual = 0 ; 
  
  if ( Steps > 0 ) {
    //
    //  Phase 1:  Compute b = A' A' A xexact
    //
    Problem.SetOperator( Amat );
   
    //
    //  We only set transpose if we have to - this allows valgrind to check
    //  that transpose is set to a default value before it is used.
    //
    if ( transpose ) OUR_CHK_ERR( Abase->SetUseTranspose( transpose ) ); 
    //    if (verbose) ParamList.set( "DebugLevel", 1 );
    //    if (verbose) ParamList.set( "OutputLevel", 1 );
    if ( Steps > 1 ) {
      OUR_CHK_ERR( Abase->SetParameters( ParamList ) ); 
      if ( Steps > 2 ) {
		
	xexact.Random();
	xexact.PutScalar(1.0);
	
	//
	//  Compute cAx = A' xexact
	//
	Amat->Multiply( transpose, xexact, b ) ;  //  b = A x2 = A A' A'' xexact

#if 0 
	std::cout << __FILE__ << "::"  << __LINE__ << "b = " << std::endl ; 
	b.Print( std::cout ) ; 
	std::cout << __FILE__ << "::"  << __LINE__ << "xexact = " << std::endl ; 
	xexact.Print( std::cout ) ; 
	std::cout << __FILE__ << "::"  << __LINE__ << "x = " << std::endl ; 
	x.Print( std::cout ) ; 
#endif
	//
	//  Phase 2:  Solve A' A' A x = b 
	//
	//
	//  Solve A sAAx = b 
	//
	Problem.SetLHS( &x );
	Problem.SetRHS( &b );
	OUR_CHK_ERR( Abase->SymbolicFactorization(  ) ); 
	if ( Steps > 2 ) {
	  OUR_CHK_ERR( Abase->SymbolicFactorization(  ) ); 
	  if ( Steps > 3 ) {
	    OUR_CHK_ERR( Abase->NumericFactorization(  ) ); 
	    if ( Steps > 4 ) {
	      OUR_CHK_ERR( Abase->NumericFactorization(  ) ); 
	      if ( Steps > 5 ) {
		OUR_CHK_ERR( Abase->Solve(  ) ); 
		if ( Steps > 6 ) {
		  OUR_CHK_ERR( Abase->Solve(  ) ); 


		  Amat->Multiply( transpose, x, bcheck ) ; //  temp = A" x2
		  
		  double norm_diff ;
		  double norm_one ;
//.........这里部分代码省略.........

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

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

  int MyPID = Comm.MyPID();

  // Number of dimension of the domain
  int space_dim = 2;
  
  // Size of each of the dimensions of the domain
  std::vector<double> brick_dim( space_dim );
  brick_dim[0] = 1.0;
  brick_dim[1] = 1.0;
  
  // Number of elements in each of the dimensions of the domain
  std::vector<int> elements( space_dim );
  elements[0] = 10;
  elements[1] = 10;
  
  // Create problem
  Teuchos::RCP<ModalProblem> testCase = Teuchos::rcp( new ModeLaplace2DQ2(Comm, brick_dim[0], elements[0], brick_dim[1], elements[1]) );
  
  // Get the stiffness and mass matrices
  Teuchos::RCP<Epetra_CrsMatrix> K = Teuchos::rcp( const_cast<Epetra_CrsMatrix *>(testCase->getStiffness()), false );
  Teuchos::RCP<Epetra_CrsMatrix> M = Teuchos::rcp( const_cast<Epetra_CrsMatrix *>(testCase->getMass()), false );
	
  //
  // *******************************************************
  // Set up Amesos direct solver for inner iteration
  // *******************************************************
  //

  // Create the shifted system K - sigma * M.
  // For the buckling transformation, this shift must be nonzero.

  double sigma = 1.0;
  Epetra_CrsMatrix Kcopy( *K );

  int addErr = EpetraExt::MatrixMatrix::Add( *M, false, -sigma, Kcopy, 1.0 );
  if (addErr != 0) {
    if (MyPID == 0) {
      std::cout << "EpetraExt::MatrixMatrix::Add returned with error: " << addErr << std::endl;
    }
#ifdef HAVE_MPI
    MPI_Finalize() ;
#endif
    return -1;
  }

  // Create Epetra linear problem class to solve "x = b"
  Epetra_LinearProblem AmesosProblem;
  AmesosProblem.SetOperator(&Kcopy);
  
  // Create Amesos factory and solver for solving "(K - sigma*M)x = b" using a direct factorization
  Amesos amesosFactory;
  Teuchos::RCP<Amesos_BaseSolver> AmesosSolver = 
    Teuchos::rcp( amesosFactory.Create( "Klu", AmesosProblem ) );

  // The AmesosBucklingOp class assumes that the symbolic/numeric factorizations have already
  // been performed on the linear problem.
  AmesosSolver->SymbolicFactorization();
  AmesosSolver->NumericFactorization();

  //
  // ************************************
  // Start the block Arnoldi iteration
  // ************************************
  //
  //  Variables used for the Block Arnoldi Method
  //
  int nev = 10;
  int blockSize = 3;  
  int numBlocks = 3*nev/blockSize;
  int maxRestarts = 5;
  //int step = 5;
  double tol = 1.0e-8;
  std::string which = "LM";
  int verbosity = Anasazi::Errors + Anasazi::Warnings + Anasazi::FinalSummary;
  //
  // 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( "Num Blocks", numBlocks );
  MyPL.set( "Maximum Restarts", maxRestarts );
  MyPL.set( "Convergence Tolerance", tol );
  //MyPL.set( "Step Size", step );
  
  typedef Epetra_MultiVector MV;
  typedef Epetra_Operator OP;
  typedef Anasazi::MultiVecTraits<double, MV> MVT;
//.........这里部分代码省略.........

double Ifpack_Condest(const Ifpack_Preconditioner& IFP,
		      const Ifpack_CondestType CT,
		      const int MaxIters,
		      const double Tol,
		      Epetra_RowMatrix* Matrix)
{
  double ConditionNumberEstimate = -1.0;

  if (CT == Ifpack_Cheap) {

    // Create a vector with all values equal to one
    Epetra_Vector Ones(IFP.OperatorDomainMap());
    Ones.PutScalar(1.0);
    // Create the vector of results
    Epetra_Vector OnesResult(IFP.OperatorRangeMap());
    // Compute the effect of the solve on the vector of ones
    IFPACK_CHK_ERR(IFP.ApplyInverse(Ones, OnesResult)); 
    // Make all values non-negative
    IFPACK_CHK_ERR(OnesResult.Abs(OnesResult)); 
    // Get the maximum value across all processors
    IFPACK_CHK_ERR(OnesResult.MaxValue(&ConditionNumberEstimate)); 

  }
  else if (CT == Ifpack_CG) {

#ifdef HAVE_IFPACK_AZTECOO
    if (Matrix == 0)
      Matrix = (Epetra_RowMatrix*)&(IFP.Matrix());

    Epetra_Vector LHS(IFP.OperatorDomainMap());
    LHS.PutScalar(0.0);
    Epetra_Vector RHS(IFP.OperatorRangeMap());
    RHS.Random();
    Epetra_LinearProblem Problem;
    Problem.SetOperator(Matrix);
    Problem.SetLHS(&LHS);
    Problem.SetRHS(&RHS);

    AztecOO Solver(Problem);
    Solver.SetAztecOption(AZ_output,AZ_none);
    Solver.SetAztecOption(AZ_solver,AZ_cg_condnum);
    Solver.Iterate(MaxIters,Tol);

    const double* status = Solver.GetAztecStatus();
    ConditionNumberEstimate = status[AZ_condnum];
#endif

  } else if (CT == Ifpack_GMRES) {

#ifdef HAVE_IFPACK_AZTECOO
    if (Matrix == 0)
      Matrix = (Epetra_RowMatrix*)&(IFP.Matrix());

    Epetra_Vector LHS(IFP.OperatorDomainMap());
    LHS.PutScalar(0.0);
    Epetra_Vector RHS(IFP.OperatorRangeMap());
    RHS.Random();
    Epetra_LinearProblem Problem;
    Problem.SetOperator(Matrix);
    Problem.SetLHS(&LHS);
    Problem.SetRHS(&RHS);

    AztecOO Solver(Problem);
    Solver.SetAztecOption(AZ_solver,AZ_gmres_condnum);
    Solver.SetAztecOption(AZ_output,AZ_none);
    // the following can be problematic for large problems,
    // but any restart would destroy useful information about
    // the condition number.
    Solver.SetAztecOption(AZ_kspace,MaxIters);
    Solver.Iterate(MaxIters,Tol);

    const double* status = Solver.GetAztecStatus();
    ConditionNumberEstimate = status[AZ_condnum];
#endif
  }

  return(ConditionNumberEstimate);

}

int
main (int argc, char *argv[])
{
  using Teuchos::RCP;
  using Teuchos::rcp;
  using std::cerr;
  using std::cout;
  using std::endl;
  // Anasazi solvers have the following template parameters:
  //
  //   - Scalar: The type of dot product results.
  //   - MV: The type of (multi)vectors.
  //   - OP: The type of operators (functions from multivector to
  //     multivector).  A matrix (like Epetra_CrsMatrix) is an example
  //     of an operator; an Ifpack preconditioner is another example.
  //
  // Here, Scalar is double, MV is Epetra_MultiVector, and OP is
  // Epetra_Operator.
  typedef Epetra_MultiVector MV;
  typedef Epetra_Operator OP;
  typedef Anasazi::MultiVecTraits<double, MV> MVT;

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

  const int MyPID = Comm.MyPID ();

  //
  // Set up the test problem
  //

  // Dimensionality of the spatial domain to discretize
  const int space_dim = 2;

  // Size of each of the dimensions of the (discrete) domain
  std::vector<double> brick_dim (space_dim);
  brick_dim[0] = 1.0;
  brick_dim[1] = 1.0;

  // Number of elements in each of the dimensions of the domain
  std::vector<int> elements (space_dim);
  elements[0] = 10;
  elements[1] = 10;

  // Create the test problem.
  RCP<ModalProblem> testCase =
    rcp (new ModeLaplace2DQ2 (Comm, brick_dim[0], elements[0],
                              brick_dim[1], elements[1]));

  // Get the stiffness and mass matrices.
  //
  // rcp (T*, false) returns a nonowning (doesn't deallocate) RCP.
  RCP<Epetra_CrsMatrix> K =
    rcp (const_cast<Epetra_CrsMatrix* > (testCase->getStiffness ()), false);
  RCP<Epetra_CrsMatrix> M =
    rcp (const_cast<Epetra_CrsMatrix* > (testCase->getMass ()), false);

  //
  // Create linear solver for linear systems with K
  //
  // Anasazi uses shift and invert, with a "shift" of zero, to find
  // the eigenvalues of least magnitude.  In this example, we
  // implement the "invert" part of shift and invert by using an
  // Amesos direct solver.
  //

  // Create Epetra linear problem class for solving linear systems
  // with K.  This implements the inverse operator for shift and
  // invert.
  Epetra_LinearProblem AmesosProblem;
  // Tell the linear problem about the matrix K.  Epetra_LinearProblem
  // doesn't know about RCP, so we have to give it a raw pointer.
  AmesosProblem.SetOperator (K.get ());

  // Create Amesos factory and solver for solving linear systems with
  // K.  The solver uses the KLU library to do a sparse LU
  // factorization.
  //
  // Note that the AmesosProblem object "absorbs" K.  Anasazi doesn't
  // see K, just the operator that implements K^{-1} M.
  Amesos amesosFactory;
  RCP<Amesos_BaseSolver> AmesosSolver;

  // Amesos can interface to many different solvers.  The following
  // loop picks a solver that Amesos supports.  The loop order
  // reflects solver preference, only in the sense that using LAPACK
  // here is a suboptimal fall-back.  (With the LAPACK option, Amesos
  // makes a dense version of the sparse matrix and uses LAPACK to
  // compute the factorization.  The other options are true sparse
  // direct factorizations.)
  const int numSolverNames = 9;
  const char* solverNames[9] = {
    "Klu", "Umfpack", "Superlu", "Superludist", "Mumps",
    "Paradiso", "Taucs", "CSparse", "Lapack"
  };
//.........这里部分代码省略.........

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

  Teuchos::GlobalMPISession mpisess(&argc,&argv,&cout);

  bool run_me = (mpisess.getRank() == 0);
  bool testFailed = false;

  if (run_me) {
    try {
      // useful typedefs
      typedef double                              ST;
      typedef Teuchos::ScalarTraits<ST>          STT;
      typedef STT::magnitudeType                  MT;
      typedef Teuchos::ScalarTraits<MT>          MTT;
      typedef Epetra_MultiVector                  MV;
      typedef Epetra_Operator                     OP;
      typedef Anasazi::MultiVecTraits<ST,MV>     MVT;
      typedef Anasazi::OperatorTraits<ST,MV,OP>  OPT;

      // parse here, so everyone knows about failure
      bool verbose    = false;
      bool debug      = false;
      std::string k_filename = "linearized_qevp_A.hb";
      std::string m_filename = "linearized_qevp_B.hb";
      int blockSize   = 1;
      int numBlocks   = 30;
      int nev         = 20;
      int maxRestarts = 0;
      int extraBlocks = 0;
      int stepSize    = 1;
      int numPrint    = 536;
      MT  tol = 1e-8;
      Teuchos::CommandLineProcessor cmdp(true,true);
      cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
      cmdp.setOption("debug","normal",&debug,"Print debugging information.");
      cmdp.setOption("A-filename",&k_filename,"Filename and path of the stiffness matrix.");
      cmdp.setOption("B-filename",&m_filename,"Filename and path of the mass matrix.");
      cmdp.setOption("extra-blocks",&extraBlocks,"Extra blocks to keep on a restart.");
      cmdp.setOption("block-size",&blockSize,"Block size.");
      cmdp.setOption("num-blocks",&numBlocks,"Number of blocks in Krylov basis.");
      cmdp.setOption("step-size",&stepSize,"Step size.");
      cmdp.setOption("nev",&nev,"Number of eigenvalues to compute.");
      cmdp.setOption("num-restarts",&maxRestarts,"Maximum number of restarts.");
      cmdp.setOption("tol",&tol,"Convergence tolerance.");
      cmdp.setOption("num-print",&numPrint,"Number of Ritz values to print.");
      // parse() will throw an exception on error
      cmdp.parse(argc,argv);

      // Get the stiffness and mass matrices
      Epetra_SerialComm Comm;
      RCP<Epetra_Map> Map;
      RCP<Epetra_CrsMatrix> A, B;
      EpetraExt::readEpetraLinearSystem( k_filename, Comm, &A, &Map );
      EpetraExt::readEpetraLinearSystem( m_filename, Comm, &B, &Map );

      //
      // *******************************************************
      // Set up Amesos direct solver for inner iteration
      // *******************************************************
      //
      // Create Epetra linear problem class to solve "Kx = b"
      Epetra_LinearProblem AmesosProblem;
      AmesosProblem.SetOperator(A.get());
      // Create Amesos factory and solver for solving "Kx = b" using a direct factorization
      Amesos amesosFactory;
      RCP<Amesos_BaseSolver> AmesosSolver = rcp( amesosFactory.Create( "Klu", AmesosProblem ) );
      // The AmesosGenOp class assumes that the symbolic/numeric factorizations have already
      // been performed on the linear problem.
      AmesosSolver->SymbolicFactorization();
      AmesosSolver->NumericFactorization();

      //
      // ************************************
      // Start the block Arnoldi iteration
      // ************************************
      //
      //  Variables used for the Block Arnoldi Method
      //
      int verbosity = Anasazi::Errors + Anasazi::Warnings + Anasazi::FinalSummary;
      if (verbose) {
        verbosity += Anasazi::IterationDetails;
      }
      if (debug) {
        verbosity += Anasazi::Debug;
      }
      //
      // Create parameter list to pass into solver
      //
      Teuchos::ParameterList MyPL;
      MyPL.set( "Verbosity", verbosity );
      MyPL.set( "Which", "LM" );
      MyPL.set( "Block Size", blockSize );
      MyPL.set( "Num Blocks", numBlocks );
      MyPL.set( "Maximum Restarts", maxRestarts );
      MyPL.set( "Convergence Tolerance", tol );
      MyPL.set( "Step Size", stepSize );
      MyPL.set( "Extra NEV Blocks", extraBlocks );
      MyPL.set( "Print Number of Ritz Values", numPrint );

      // Create an Epetra_MultiVector for an initial vector to start the solver.
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
#ifdef HAVE_MPI
    Teuchos::GlobalMPISession mpiSession(&argc, &argv, 0);
    Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
    Epetra_SerialComm Comm;
#endif
    int nProcs, myPID ;
    Teuchos::ParameterList pLUList ;        // ParaLU parameters
    Teuchos::ParameterList isoList ;        // Isorropia parameters
    string ipFileName = "ShyLU.xml";       // TODO : Accept as i/p

    nProcs = mpiSession.getNProc();
    myPID = Comm.MyPID();

    if (myPID == 0)
    {
        cout <<"Parallel execution: nProcs="<< nProcs << endl;
    }

    // =================== Read input xml file =============================
    Teuchos::updateParametersFromXmlFile(ipFileName, &pLUList);
    isoList = pLUList.sublist("Isorropia Input");
    // Get matrix market file name
    string MMFileName = Teuchos::getParameter<string>(pLUList, "mm_file");
    string prec_type = Teuchos::getParameter<string>(pLUList, "preconditioner");

    if (myPID == 0)
    {
        cout << "Input :" << endl;
        cout << "ParaLU params " << endl;
        pLUList.print(std::cout, 2, true, true);
        cout << "Matrix market file name: " << MMFileName << endl;
    }

    // ==================== Read input Matrix ==============================
    Epetra_CrsMatrix *A;

    int err = EpetraExt::MatrixMarketFileToCrsMatrix(MMFileName.c_str(), Comm, A);
    //EpetraExt::MatlabFileToCrsMatrix(MMFileName.c_str(), Comm, A);
    //assert(err != 0);
    cout <<"Done reading the matrix"<< endl;
    int n = A->NumGlobalRows();
    cout <<"n="<< n << endl;

    // Create input vectors
    Epetra_Map vecMap(n, 0, Comm);
    Epetra_MultiVector x(vecMap, 1);
    Epetra_MultiVector b(vecMap, 1, false);
    b.PutScalar(1.0); // TODO : Accept it as input

    // Partition the matrix with hypergraph partitioning and redisstribute
    Isorropia::Epetra::Partitioner *partitioner = new
                            Isorropia::Epetra::Partitioner(A, isoList, false);
    partitioner->partition();
    Isorropia::Epetra::Redistributor rd(partitioner);

    Epetra_CrsMatrix *newA;
    Epetra_MultiVector *newX, *newB; 
    rd.redistribute(*A, newA);
    delete A;
    A = newA;

    rd.redistribute(x, newX);
    rd.redistribute(b, newB);

    Epetra_LinearProblem problem(A, newX, newB);

    Amesos Factory;
    char* SolverType = "Amesos_Klu";
    bool IsAvailable = Factory.Query(SolverType);

    Epetra_LinearProblem *LP = new Epetra_LinearProblem();
    LP->SetOperator(A);
    LP->SetLHS(newX);
    LP->SetRHS(newB);
    Amesos_BaseSolver *Solver = Factory.Create(SolverType, *LP);


    Solver->SymbolicFactorization();
  Teuchos::Time ftime("setup time");
      ftime.start();
    Solver->NumericFactorization();
    cout << "Numeric Factorization" << endl;
    Solver->Solve();
    cout << "Solve done" << endl;

    ftime.stop();
    cout << "Time to setup" << ftime.totalElapsedTime() << endl;

    // compute ||Ax - b||
    double Norm;
    Epetra_MultiVector Ax(vecMap, 1);

    Epetra_MultiVector *newAx; 
    rd.redistribute(Ax, newAx);
    A->Multiply(false, *newX, *newAx);
    newAx->Update(1.0, *newB, -1.0);
    newAx->Norm2(&Norm);
//.........这里部分代码省略.........