C++ Amesos::Create方法代码示例

  void AmesosSmoother<Node>::Setup(Level& currentLevel) {
    FactoryMonitor m(*this, "Setup Smoother", currentLevel);

    if (SmootherPrototype::IsSetup() == true)
      this->GetOStream(Warnings0) << "MueLu::AmesosSmoother::Setup(): Setup() has already been called" << std::endl;

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

    RCP<Epetra_CrsMatrix> epA = Utilities::Op2NonConstEpetraCrs(A_);
    linearProblem_ = rcp( new Epetra_LinearProblem() );
    linearProblem_->SetOperator(epA.get());

    Amesos factory;
    prec_ = rcp(factory.Create(type_, *linearProblem_));
    TEUCHOS_TEST_FOR_EXCEPTION(prec_ == Teuchos::null, Exceptions::RuntimeError, "MueLu::AmesosSmoother::Setup(): Solver '" + type_ + "' not supported by Amesos");

    // set Reindex flag, if A is distributed with non-contiguous maps
    // unfortunately there is no reindex for Amesos2, yet. So, this only works for Epetra based problems
    if (A_->getRowMap()->isDistributed() == true && A_->getRowMap()->isContiguous() == false)
      const_cast<ParameterList&>(this->GetParameterList()).set("Reindex", true);

    const ParameterList& paramList = this->GetParameterList();
    RCP<ParameterList> precList = this->RemoveFactoriesFromList(paramList);

    prec_->SetParameters(*precList);

    const_cast<ParameterList&>(paramList).setParameters(*precList);

    int r = prec_->NumericFactorization();
    TEUCHOS_TEST_FOR_EXCEPTION(r != 0, Exceptions::RuntimeError, "MueLu::AmesosSmoother::Setup(): Amesos solver returns value of " +
                               Teuchos::toString(r) + " during NumericFactorization()");

    SmootherPrototype::IsSetup(true);
  }

void DistributedOverlapMatrix<Rank>::SolveOverlapRank(Wavefunction<Rank> &srcPsi, Wavefunction<Rank> &destPsi, int opRank, bool fastAlgo)
//void DistributedOverlapMatrix<Rank>::SolveOverlapRank(Wavefunction<Rank> &psi, int opRank, bool fastAlgo)
{
	SetupRank(srcPsi, opRank);
	
	if (fastAlgo)
	{
		//Setup source multivector
		WavefunctionToMultiVector(srcPsi, InputVector(opRank), opRank);
		
		//Solve
		AMESOS_CHK_ERRV(Solvers(opRank)->Solve());

		//Put solution multivector into destination wavefunction
		MultiVectorToWavefunction(destPsi, OutputVector(opRank), opRank);
	}
	else
	{
		//Setup source multivector
		Epetra_MultiVector_Ptr srcVec =  SetupMultivector(srcPsi, opRank);

		//Output multivector (copy of input)
		Epetra_MultiVector_Ptr destVec = Epetra_MultiVector_Ptr(new Epetra_MultiVector(*srcVec));

		//Set up Epetra LinearProblem with overlap for this rank and input/output multivectors
		Epetra_LinearProblem Problem(OverlapMatrices(opRank).get(), destVec.get(), srcVec.get());

		//Initialize Amesos solver and factory
		Amesos_BaseSolver* Solver;
		Amesos Factory;
		std::string SolverType = "Amesos_Superludist";
		Solver = Factory.Create(SolverType, Problem);

		//Check that requested solver exists in Amesos
		if (Solver == 0)
		{
			throw std::runtime_error("Specified Amesos solver not available");
		}

		//Setup the parameter list for the solver
		Teuchos::ParameterList List;
		//List.set("PrintTiming", true);
		//List.set("PrintStatus", true);
		Solver->SetParameters(List);

		//Factorization and solve
		Solver->SymbolicFactorization();
		Solver->NumericFactorization();
		AMESOS_CHK_ERRV(Solver->Solve());

		//Put solution multivector into destination wavefunction
		MultiVectorToWavefunction(destPsi, destVec, opRank);

		//Solver->PrintTiming();
	}

}

SolverState<double> AmesosSolver::solve(const LinearOperator<double>& op, 
  const Vector<double>& rhs, 
  Vector<double>& soln) const
{
	Playa::Vector<double> bCopy = rhs.copy();
	Playa::Vector<double> xCopy = rhs.copy();

  Epetra_Vector* b = EpetraVector::getConcretePtr(bCopy);
  Epetra_Vector* x = EpetraVector::getConcretePtr(xCopy);

	Epetra_CrsMatrix& A = EpetraMatrix::getConcrete(op);

  Epetra_LinearProblem prob(&A, x, b);

  Amesos amFactory;
  RCP<Amesos_BaseSolver> solver 
    = rcp(amFactory.Create("Amesos_" + kernel_, prob));
  TEUCHOS_TEST_FOR_EXCEPTION(solver.get()==0, std::runtime_error, 
    "AmesosSolver::solve() failed to instantiate "
    << kernel_ << "solver kernel");

  int ierr = solver->Solve();
  
  soln = xCopy;

  SolverStatusCode state;
  std::string msg;

  switch(ierr)
  {
    case 0:
      state = SolveConverged;
      msg = "converged";
      break;
    default:
      state = SolveCrashed;
      msg = "amesos failed: ierr=" + Teuchos::toString(ierr);
  }

  SolverState<double> rtn(state, "Amesos solver " + msg, 0, 0);
  return rtn;
}

void DistributedOverlapMatrix<Rank>::SetupOverlapSolvers(Wavefunction<Rank> &psi, int opRank)
{
	//Set up Epetra LinearProblem with overlap for this rank and input/output multivectors
	EpetraProblems(opRank) = Epetra_LinearProblem_Ptr( new Epetra_LinearProblem(OverlapMatrices(opRank).get(), OutputVector(opRank).get(), InputVector(opRank).get()) );

	//Determine solver type. Use SuperLU_dist if opRank is distributed, else KLU
	Amesos Factory;
	std::string SolverType;
	if (psi.GetRepresentation()->GetDistributedModel()->IsDistributedRank(opRank))
	{
		SolverType = "Amesos_Superludist";
	}
	else
	{
		SolverType = "Amesos_Klu";
	}

	//Create Amesos solver for this rank
	Solvers(opRank) = Amesos_BaseSolver_Ptr( Factory.Create(SolverType, *EpetraProblems(opRank)) );

	//Check that requested solver exists in Amesos build
	if (Solvers(opRank) == 0)
	{
		throw std::runtime_error("Specified Amesos solver not available");
	}

	//Setup the parameter list for the solver (TODO: get these from Python config section)
	Teuchos::ParameterList List;
	List.set("MatrixType", "symmetric");
	List.set("PrintTiming", false);
	List.set("PrintStatus", false);
	List.set("ComputeTrueResidual", false);
	Solvers(opRank)->SetParameters(List);

	//Perform symbolic factorization
	Solvers(opRank)->SymbolicFactorization();
	Solvers(opRank)->NumericFactorization();
}

//-----------------------------------------------------------------------------
// Function      : N_LAS_HBBlockJacobiPrecond::initGraph
// Purpose       : Initialize the graph using information from the N_LAS_System
// Special Notes :
// Scope         : Public
// Creator       : Heidi Thornquist, SNL, Electrical & Microsystem Modeling
// Creation Date : 12/11/08
//-----------------------------------------------------------------------------
bool N_LAS_HBBlockJacobiPrecond::initGraph( const Teuchos::RCP<N_LAS_Problem> & problem )
{
  // Generate the graph of each real equivalent form and then generate
  // empty linear systems for each frequency.
  RCP<N_LAS_Matrix> appdQdx = rcp( appBuilderPtr_->createMatrix() );
  RCP<N_LAS_Matrix> appdFdx = rcp( appBuilderPtr_->createMatrix() );

  // Relate the conductance and inductance matrices:
  // G(t) = df/dx(x(t))
  // C(t) = dq/dx(x(t))
  // Compute:  (omega*i*j*C_bar + G_bar)^{-1} for i=0,1,...,M,-M,...,-1
  //           using real equivalent form K_1 = [G_bar -i*omega*C_bar; i*omega*C_bar G_bar]

  // Generate the new real equivalent graph
  RCP<const Epetra_Map> origMap = appBuilderPtr_->getSolutionMap(); 
  int origLocalRows = origMap->NumMyElements();
  int origGlobalRows = origMap->NumGlobalElements();
  int refRows = 2*origLocalRows;
  std::vector<int> rowIdxs( refRows );
  int * origIdxs = origMap->MyGlobalElements();
  for (int i=0; i<origLocalRows; ++i)
  {
    rowIdxs[i] = origIdxs[i];
    rowIdxs[origLocalRows+i] = origIdxs[i] + origGlobalRows;
  }
  epetraMap_ = rcp(new Epetra_Map( -1, refRows, &rowIdxs[0], 0, (appdQdx->epetraObj()).Comm() ) );

  // Count up the number of nonzero entries for the 2x2 block matrix.
  std::vector<int> refNNZs(refRows);
  maxRefNNZs_ = 0;
  for ( int i=0; i<origLocalRows; ++i ) {
    refNNZs[i] = appdQdx->getLocalRowLength(i) + appdFdx->getLocalRowLength(i);
    refNNZs[origLocalRows+i] = refNNZs[i];
    if (refNNZs[i] > maxRefNNZs_) maxRefNNZs_ = refNNZs[i];
  }
  epetraGraph_ = rcp(new Epetra_CrsGraph( Copy, *epetraMap_, &refNNZs[0], true ));

  // Put together the indices for each row and insert them into the graph.
  int tmpNNZs=0, tmpNNZs2=0;
  std::vector<double> tmpCoeffs(maxRefNNZs_);
  std::vector<int> refIdxs(maxRefNNZs_), refIdxs2(maxRefNNZs_);

  for ( int i=0; i<origLocalRows; ++i ) {

    // Get the indices for the first block of the matrix (G_bar)
    appdFdx->getRowCopy( rowIdxs[i], maxRefNNZs_, tmpNNZs, &tmpCoeffs[0], &refIdxs[0] );

    // Get the indices for the third block of the matrix (C_bar)
    appdQdx->getRowCopy( rowIdxs[i], maxRefNNZs_, tmpNNZs2, &tmpCoeffs[0], &refIdxs2[0] );

    // Insert the indices for the third block into refIdxs, as they are the indices of the second block
    for (int j=0; j<tmpNNZs2; ++j) {
      refIdxs[tmpNNZs+j] = refIdxs2[j]+origGlobalRows;
    }
    epetraGraph_->InsertGlobalIndices( rowIdxs[i], refNNZs[i], &refIdxs[0] );

    // Insert the indices for the first block into refIdxs2, as they are the indices of the fourth block
    for (int j=0; j<tmpNNZs; ++j) {
      refIdxs2[tmpNNZs2+j] = refIdxs[j]+origGlobalRows;
    }
    epetraGraph_->InsertGlobalIndices( rowIdxs[origLocalRows+i], refNNZs[origLocalRows+i], &refIdxs2[0] );
  }
  epetraGraph_->FillComplete();

  // Get the Fourier series information and generate the Epetra_LinearSystems.
  RCP<N_LAS_BlockVector> bXt = hbBuilderPtr_->createTimeDomainBlockVector();
  N_ = bXt->blockCount();
  M_ = (int)((N_-1)/2);

  // Generate the vectors for the N_ linear problems to be solved on the diagonal.
  epetraRHS_ = rcp( new Epetra_MultiVector( *epetraMap_, 1 ) );
  epetraSoln_ = rcp( new Epetra_MultiVector( *epetraMap_, 1 ) );
  epetraMatrix_.resize(N_);
  epetraProblem_.resize(N_);
  amesosPtr_.resize(N_);

  Amesos amesosFactory;
  Teuchos::ParameterList params;
  
#ifndef Xyce_PARALLEL_MPI 
  // Inform solver not to check inputs to reduce overhead.
  params.set( "TrustMe", true );
#endif

  for (int i=0; i<N_; ++i) {
    epetraMatrix_[i] = rcp( new Epetra_CrsMatrix( Copy, *epetraGraph_ ) );
    epetraMatrix_[i]->FillComplete();
    epetraMatrix_[i]->OptimizeStorage();
    epetraProblem_[i] = rcp( new Epetra_LinearProblem( &*epetraMatrix_[i], &*epetraRHS_, &*epetraSoln_ ) );
    amesosPtr_[i] = rcp( amesosFactory.Create( "Klu", *epetraProblem_[i] ) );
    amesosPtr_[i]->SetParameters( params );
    amesosPtr_[i]->SymbolicFactorization();
//.........这里部分代码省略.........

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

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

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

  bool verbose = (Comm.MyPID() == 0);
  double TotalResidual = 0.0;

  // Create the Map, defined as a grid, of size nx x ny x nz,
  // subdivided into mx x my x mz cubes, each assigned to a 
  // different processor.

#if 0
  ParameterList GaleriList;
  GaleriList.set("nx", 4);
  GaleriList.set("ny", 4);
  GaleriList.set("nz", 4 * Comm.NumProc());
  GaleriList.set("mx", 1);
  GaleriList.set("my", 1);
  GaleriList.set("mz", Comm.NumProc());
  Epetra_Map* Map = CreateMap("Cartesian3D", Comm, GaleriList);
  
  // Create a matrix, in this case corresponding to a 3D Laplacian
  // discretized using a classical 7-point stencil. Please refer to 
  // the Galeri documentation for an overview of available matrices.
  // 
  // NOTE: matrix must be symmetric if DSCPACK is used.

  Epetra_CrsMatrix* Matrix = CreateCrsMatrix("Laplace3D", Map, GaleriList);
#else
  bool transpose = false ; 
  bool distribute = false ; 
  bool symmetric ; 
  Epetra_CrsMatrix *Matrix = 0 ;
  Epetra_Map *Map = 0 ;
  CreateCrsMatrix( "ibm.triU", Comm, Map, transpose, distribute, &symmetric, Matrix ) ;




#endif

  // build vectors, in this case with 1 vector
  Epetra_MultiVector LHS(*Map, 1); 
  Epetra_MultiVector RHS(*Map, 1); 
    
  // create a linear problem object
  Epetra_LinearProblem Problem(Matrix, &LHS, &RHS);

  // use this list to set up parameters, now it is required
  // to use all the available processes (if supported by the
  // underlying solver). Uncomment the following two lines
  // to let Amesos print out some timing and status information.
  ParameterList List;
  List.set("PrintTiming",true);
  List.set("PrintStatus",true);
  List.set("MaxProcs",Comm.NumProc());

  std::vector<std::string> SolverType;
  SolverType.push_back("Amesos_Paraklete");

  Epetra_Time Time(Comm);
  
  // this is the Amesos factory object that will create 
  // a specific Amesos solver.
  Amesos Factory;

  // Cycle over all solvers.
  // Only installed solvers will be tested.
  for (unsigned int i = 0 ; i < SolverType.size() ; ++i) 
  {
    // Check whether the solver is available or not
    if (Factory.Query(SolverType[i])) 
    {
      // 1.- set exact solution (constant vector)
      LHS.PutScalar(1.0);
 
      // 2.- create corresponding rhs
      Matrix->Multiply(false, LHS, RHS);
 
      // 3.- randomize solution vector
      LHS.Random();
 
      // 4.- create the amesos solver object
      Amesos_BaseSolver* Solver = Factory.Create(SolverType[i], Problem);
      assert (Solver != 0);

      Solver->SetParameters(List);

      // 5.- factorize and solve
      
      Time.ResetStartTime();
      AMESOS_CHK_ERR(Solver->SymbolicFactorization());

      // 7.- delete the object
      delete Solver;
//.........这里部分代码省略.........

//==============================================================================
int Ifpack2_Amesos::Initialize()
{

  IsInitialized_ = false;
  IsComputed_ = false;

  if (Matrix_ == Teuchos::null)
    IFPACK2_CHK_ERR(-1);

#if 0
  // better to avoid strange games with maps, this class should be
  // used for Ifpack2_LocalFilter'd matrices only
  if (Comm().NumProc() != 1) {
    cout << "Class Ifpack2_Amesos must be used for serial runs;" << endl;
    cout << "for parallel runs you should declare objects as:" << endl; 
    cout << "Ifpack2_AdditiveSchwarz<Ifpack2_Amesos> APrec(Matrix)" << endl;
    exit(EXIT_FAILURE);
  }
#endif

  // only square matrices
  if (Matrix_->NumGlobalRows() != Matrix_->NumGlobalCols())
    IFPACK2_CHK_ERR(-1);

  // at least one nonzero
  if (Matrix_->NumMyNonzeros() == 0) 
    IFPACK2_CHK_ERR(-1);

  Problem_->SetOperator(const_cast<Tpetra_RowMatrix*>(Matrix_.get()));

  if (Time_ == Teuchos::null)
    Time_ = Teuchos::rcp( new Tpetra_Time(Comm()) );

  Amesos Factory;
  Solver_ = Teuchos::rcp( Factory.Create((char*)Label_.c_str(),*Problem_) );
  
  if (Solver_ == Teuchos::null) 
  {
    // try to create KLU, it is generally enabled
    Solver_ = Teuchos::rcp( Factory.Create("Amesos_Klu",*Problem_) );
  }
  if (Solver_ == Teuchos::null)
  {
    // finally try to create LAPACK, it is generally enabled
    // NOTE: the matrix is dense, so this should only be for
    // small problems...
    if (FirstTime)
    {
      cerr << "TIFPACK WARNING: In class Ifpack2_Amesos:" << endl;
      cerr << "TIFPACK WARNING: Using LAPACK because other Amesos" << endl;
      cerr << "TIFPACK WARNING: solvers are not available. LAPACK" << endl;
      cerr << "TIFPACK WARNING: allocates memory to store the matrix as" << endl;
      cerr << "TIFPACK WARNING: dense, I hope you have enough memory..." << endl;
      cerr << "TIFPACK WARNING: (file " << __FILE__ << ", line " << __LINE__ 
           << ")" << endl;
      FirstTime = false;
    }
    Solver_ = Teuchos::rcp( Factory.Create("Amesos_Lapack",*Problem_) );
  }
  // if empty, give up.
  if (Solver_ == Teuchos::null)
    IFPACK2_CHK_ERR(-1);

  IFPACK2_CHK_ERR(Solver_->SetUseTranspose(UseTranspose_));
  Solver_->SetParameters(List_);
  IFPACK2_CHK_ERR(Solver_->SymbolicFactorization());

  IsInitialized_ = true;
  ++NumInitialize_;
  InitializeTime_ += Time_->ElapsedTime();
  return(0);
}

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

  // 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"
  };
  for (int k = 0; k < numSolverNames; ++k) {
    const char* const solverName = solverNames[k];
    if (amesosFactory.Query (solverName)) {
      AmesosSolver = rcp (amesosFactory.Create (solverName, AmesosProblem));
      if (MyPID == 0) {
        cout << "Amesos solver: \"" << solverName << "\"" << endl;
      }
    }
  }
  if (AmesosSolver.is_null ()) {
    throw std::runtime_error ("Amesos appears not to have any solvers enabled.");
  }

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

  //
  // Set parameters for the block Krylov-Schur eigensolver
  //

  double tol = 1.0e-8; // convergence tolerance
  int nev = 10; // number of eigenvalues for which to solve
  int blockSize = 3; // block size (number of eigenvectors processed at once)
  int numBlocks = 3 * nev / blockSize; // restart length
  int maxRestarts = 5; // maximum number of restart cycles

  // We're looking for the largest-magnitude eigenvalues of the
  // _inverse_ operator, thus, the smallest-magnitude eigenvalues of
  // the original operator.
  std::string which = "LM";
  int verbosity = Anasazi::Errors + Anasazi::Warnings + Anasazi::FinalSummary;

  // Create ParameterList to pass into eigensolver
  Teuchos::ParameterList MyPL;

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

  bool verbose = (Comm.MyPID() == 0);
  double TotalResidual = 0.0;

  // Create the Map, defined as a grid, of size nx x ny x nz,
  // subdivided into mx x my x mz cubes, each assigned to a
  // different processor.

#ifndef FILENAME_SPECIFIED_ON_COMMAND_LINE
  ParameterList GaleriList;
  GaleriList.set("nx", 4);
  GaleriList.set("ny", 4);
  GaleriList.set("nz", 4 * Comm.NumProc());
  GaleriList.set("mx", 1);
  GaleriList.set("my", 1);
  GaleriList.set("mz", Comm.NumProc());
  Epetra_Map* Map = CreateMap("Cartesian3D", Comm, GaleriList);

  // Create a matrix, in this case corresponding to a 3D Laplacian
  // discretized using a classical 7-point stencil. Please refer to
  // the Galeri documentation for an overview of available matrices.
  //
  // NOTE: matrix must be symmetric if DSCPACK is used.

  Epetra_CrsMatrix* Matrix = CreateCrsMatrix("Laplace3D", Map, GaleriList);
#else
  bool transpose = false ;
  bool distribute = false ;
  bool symmetric ;
  Epetra_CrsMatrix *Matrix = 0 ;
  Epetra_Map *Map = 0 ;
  MyCreateCrsMatrix( argv[1], Comm, Map, transpose, distribute, symmetric, Matrix ) ;

#endif

  // build vectors, in this case with 1 vector
  Epetra_MultiVector LHS(*Map, 1);
  Epetra_MultiVector RHS(*Map, 1);

  // create a linear problem object
  Epetra_LinearProblem Problem(Matrix, &LHS, &RHS);

  // use this list to set up parameters, now it is required
  // to use all the available processes (if supported by the
  // underlying solver). Uncomment the following two lines
  // to let Amesos print out some timing and status information.
  ParameterList List;
  List.set("PrintTiming",true);
  List.set("PrintStatus",true);
  List.set("MaxProcs",Comm.NumProc());

  std::vector<std::string> SolverType;
  SolverType.push_back("Amesos_Paraklete");
  SolverType.push_back("Amesos_Klu");
  Comm.Barrier() ;
#if 1
  SolverType.push_back("Amesos_Lapack");
  SolverType.push_back("Amesos_Umfpack");
  SolverType.push_back("Amesos_Pardiso");
  SolverType.push_back("Amesos_Taucs");
  SolverType.push_back("Amesos_Superlu");
  SolverType.push_back("Amesos_Superludist");
  SolverType.push_back("Amesos_Mumps");
  SolverType.push_back("Amesos_Dscpack");
  SolverType.push_back("Amesos_Scalapack");
#endif
  Epetra_Time Time(Comm);

  // this is the Amesos factory object that will create
  // a specific Amesos solver.
  Amesos Factory;

  // Cycle over all solvers.
  // Only installed solvers will be tested.
  for (unsigned int i = 0 ; i < SolverType.size() ; ++i)
  {
    // Check whether the solver is available or not
    if (Factory.Query(SolverType[i]))
    {
      // 1.- set exact solution (constant vector)
      LHS.PutScalar(1.0);

      // 2.- create corresponding rhs
      Matrix->Multiply(false, LHS, RHS);

      // 3.- randomize solution vector
      LHS.Random();

      // 4.- create the amesos solver object
      Amesos_BaseSolver* Solver = Factory.Create(SolverType[i], Problem);
      assert (Solver != 0);

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

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[]) {
#if defined(HAVE_MUELU_EPETRA)
#include <MueLu_UseShortNames.hpp>

  using Teuchos::RCP; // reference count pointers
  using Teuchos::rcp;
  using Teuchos::TimeMonitor;

  // =========================================================================
  // MPI initialization using Teuchos
  // =========================================================================
  Teuchos::GlobalMPISession mpiSession(&argc, &argv, NULL);

  bool success = false;
  try {
    RCP< const Teuchos::Comm<int> > comm = Teuchos::DefaultComm<int>::getComm();
    int MyPID   = comm->getRank();
    int NumProc = comm->getSize();

    const Teuchos::RCP<Epetra_Comm> epComm = Teuchos::rcp_const_cast<Epetra_Comm>(Xpetra::toEpetra(comm));

    // =========================================================================
    // Convenient definitions
    // =========================================================================
    //SC zero = Teuchos::ScalarTraits<SC>::zero(), one = Teuchos::ScalarTraits<SC>::one();

    // Instead of checking each time for rank, create a rank 0 stream
    RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
    Teuchos::FancyOStream& fancyout = *fancy;
    fancyout.setOutputToRootOnly(0);

    // =========================================================================
    // Parameters initialization
    // =========================================================================
    Teuchos::CommandLineProcessor clp(false);
    GO nx = 100, ny = 100;
    GO maxCoarseSize = 10;
    LO maxLevels = 4;
    clp.setOption("nx",                   &nx,              "mesh size in x direction");
    clp.setOption("ny",                   &ny,              "mesh size in y direction");
    clp.setOption("maxCoarseSize",        &maxCoarseSize,   "maximum coarse size");
    clp.setOption("maxLevels",            &maxLevels,       "maximum number of multigrid levels");
    int mgridSweeps = 1; clp.setOption("mgridSweeps", &mgridSweeps, "number of multigrid sweeps within Multigrid solver.");
    std::string printTimings = "no";   clp.setOption("timings", &printTimings,     "print timings to screen [yes/no]");
    double tol               = 1e-12;  clp.setOption("tol",                   &tol,              "solver convergence tolerance");

    switch (clp.parse(argc,argv)) {
      case Teuchos::CommandLineProcessor::PARSE_HELP_PRINTED:        return EXIT_SUCCESS; break;
      case Teuchos::CommandLineProcessor::PARSE_ERROR:
      case Teuchos::CommandLineProcessor::PARSE_UNRECOGNIZED_OPTION: return EXIT_FAILURE; break;
      case Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL:                               break;
    }

    // =========================================================================
    // Problem construction
    // =========================================================================
    RCP<TimeMonitor> globalTimeMonitor = rcp(new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: S - Global Time"))), tm;

    comm->barrier();
    tm = rcp(new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: 1 - Matrix Build")));

    Teuchos::ParameterList GaleriList;
    GaleriList.set("nx", nx);
    GaleriList.set("ny", ny);
    GaleriList.set("mx", epComm->NumProc());
    GaleriList.set("my", 1);
    GaleriList.set("lx", 1.0); // length of x-axis
    GaleriList.set("ly", 1.0); // length of y-axis
    GaleriList.set("diff", 1e-5);
    GaleriList.set("conv", 1.0);

    // create map
    Teuchos::RCP<Epetra_Map> epMap = Teuchos::rcp(Galeri::CreateMap("Cartesian2D", *epComm, GaleriList));

    // create coordinates
    Teuchos::RCP<Epetra_MultiVector> epCoord = Teuchos::rcp(Galeri::CreateCartesianCoordinates("2D", epMap.get(), GaleriList));

    // create matrix
    Teuchos::RCP<Epetra_CrsMatrix> epA = Teuchos::rcp(Galeri::CreateCrsMatrix("Recirc2D", epMap.get(), GaleriList));

    // Epetra -> Xpetra
    Teuchos::RCP<CrsMatrix> exA = Teuchos::rcp(new Xpetra::EpetraCrsMatrixT<int,Node>(epA));
    Teuchos::RCP<CrsMatrixWrap> exAWrap = Teuchos::rcp(new CrsMatrixWrap(exA));

    RCP<Matrix> A = Teuchos::rcp_dynamic_cast<Matrix>(exAWrap);
    int numPDEs = 1;
    A->SetFixedBlockSize(numPDEs);

    // set rhs and solution vector
    RCP<Epetra_Vector> B = Teuchos::rcp(new Epetra_Vector(*epMap));
    RCP<Epetra_Vector> X = Teuchos::rcp(new Epetra_Vector(*epMap));
    B->PutScalar(1.0);
    X->PutScalar(0.0);

    // Epetra -> Xpetra
    RCP<Vector> xB = Teuchos::rcp(new Xpetra::EpetraVectorT<int,Node>(B));
    RCP<Vector> xX = Teuchos::rcp(new Xpetra::EpetraVectorT<int,Node>(X));

    xX->setSeed(100);
    xX->randomize();
//.........这里部分代码省略.........

int main(int argc, char *argv[]) {
#include <MueLu_UseShortNames.hpp>

  using Teuchos::RCP; // reference count pointers
  using Teuchos::rcp;
  using Teuchos::TimeMonitor;

  // =========================================================================
  // MPI initialization using Teuchos
  // =========================================================================
  Teuchos::GlobalMPISession mpiSession(&argc, &argv, NULL);
  RCP< const Teuchos::Comm<int> > comm = Teuchos::DefaultComm<int>::getComm();
  int MyPID   = comm->getRank();
  int NumProc = comm->getSize();

  const Teuchos::RCP<Epetra_Comm> epComm = Teuchos::rcp_const_cast<Epetra_Comm>(Xpetra::toEpetra(comm));

  // =========================================================================
  // Convenient definitions
  // =========================================================================
  //SC zero = Teuchos::ScalarTraits<SC>::zero(), one = Teuchos::ScalarTraits<SC>::one();

  // Instead of checking each time for rank, create a rank 0 stream
  RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
  Teuchos::FancyOStream& fancyout = *fancy;
  fancyout.setOutputToRootOnly(0);



  // =========================================================================
  // Parameters initialization
  // =========================================================================
  Teuchos::CommandLineProcessor clp(false);
  GO nx = 100, ny = 100;
  clp.setOption("nx",                   &nx,              "mesh size in x direction");
  clp.setOption("ny",                   &ny,              "mesh size in y direction");
  std::string xmlFileName = "xml/s3a.xml"; clp.setOption("xml", &xmlFileName,     "read parameters from a file. Otherwise, this example uses by default 'tutorial1a.xml'");
  int mgridSweeps = 1; clp.setOption("mgridSweeps", &mgridSweeps, "number of multigrid sweeps within Multigrid solver.");
  std::string printTimings = "no";   clp.setOption("timings", &printTimings,     "print timings to screen [yes/no]");
  double tol               = 1e-12;  clp.setOption("tol",                   &tol,              "solver convergence tolerance");

  switch (clp.parse(argc,argv)) {
    case Teuchos::CommandLineProcessor::PARSE_HELP_PRINTED:        return EXIT_SUCCESS; break;
    case Teuchos::CommandLineProcessor::PARSE_ERROR:
    case Teuchos::CommandLineProcessor::PARSE_UNRECOGNIZED_OPTION: return EXIT_FAILURE; break;
    case Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL:                               break;
  }

  // =========================================================================
  // Problem construction
  // =========================================================================
  RCP<TimeMonitor> globalTimeMonitor = rcp(new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: S - Global Time"))), tm;

  comm->barrier();
  tm = rcp(new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: 1 - Matrix Build")));

  Teuchos::ParameterList GaleriList;
  GaleriList.set("nx", nx);
  GaleriList.set("ny", ny);
  GaleriList.set("mx", epComm->NumProc());
  GaleriList.set("my", 1);
  GaleriList.set("lx", 1.0); // length of x-axis
  GaleriList.set("ly", 1.0); // length of y-axis
  GaleriList.set("diff", 1e-5);
  GaleriList.set("conv", 1.0);

  // create map
  Teuchos::RCP<Epetra_Map> epMap = Teuchos::rcp(Galeri::CreateMap("Cartesian2D", *epComm, GaleriList));

  // create coordinates
  Teuchos::RCP<Epetra_MultiVector> epCoord = Teuchos::rcp(Galeri::CreateCartesianCoordinates("2D", epMap.get(), GaleriList));

  // create matrix
  Teuchos::RCP<Epetra_CrsMatrix> epA = Teuchos::rcp(Galeri::CreateCrsMatrix("Recirc2D", epMap.get(), GaleriList));

  // Epetra -> Xpetra
  Teuchos::RCP<CrsMatrix> exA = Teuchos::rcp(new Xpetra::EpetraCrsMatrix(epA));
  Teuchos::RCP<CrsMatrixWrap> exAWrap = Teuchos::rcp(new CrsMatrixWrap(exA));

  RCP<Matrix> A = Teuchos::rcp_dynamic_cast<Matrix>(exAWrap);
  int numPDEs = 1;
  A->SetFixedBlockSize(numPDEs);

  // set rhs and solution vector
  RCP<Epetra_Vector> B = Teuchos::rcp(new Epetra_Vector(*epMap));
  RCP<Epetra_Vector> X = Teuchos::rcp(new Epetra_Vector(*epMap));
  B->PutScalar(1.0);
  X->PutScalar(0.0);

  // Epetra -> Xpetra
  RCP<Vector> xB = Teuchos::rcp(new Xpetra::EpetraVector(B));
  RCP<Vector> xX = Teuchos::rcp(new Xpetra::EpetraVector(X));

  xX->setSeed(100);
  xX->randomize();

  // build null space vector
  RCP<const Map> map = A->getRowMap();
  RCP<MultiVector> nullspace = MultiVectorFactory::Build(map, numPDEs);

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