C++ Amesos类代码示例

  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);
  }

int ML_isKLUAvailable()
{
  Amesos Factory;
  if (Factory.Query("Amesos_Klu"))
    return(1); // this is true
  else
    return(0); // this is false
}

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();
}

  Implementation(common::Component& self) :
    m_self(self),
    m_ml_parameter_list(Teuchos::createParameterList()),
    m_solver_parameter_list(Teuchos::createParameterList()),
    m_xcoords(0),
    m_dim(0)
  {
    // ML default parameters
    ML_Epetra::SetDefaults("SA", *m_ml_parameter_list);
    m_ml_parameter_list->set("ML output", 10);
    m_ml_parameter_list->set("max levels",3);
    m_ml_parameter_list->set("aggregation: type", "METIS");
    m_ml_parameter_list->set("aggregation: nodes per aggregate", 16);
    m_ml_parameter_list->set("smoother: type","Chebyshev");
    m_ml_parameter_list->set("smoother: sweeps",2);
    m_ml_parameter_list->set("smoother: pre or post", "both");
    
    Amesos amesos;
    if(amesos.Query("Amesos_Mumps"))
    {
      m_ml_parameter_list->set("coarse: type","Amesos-MUMPS");
    }
    else
    {
      m_ml_parameter_list->set("coarse: type","Amesos-KLU");
    }
    
    // Default solver parameters
    m_solver_parameter_list->set( "Verbosity", Belos::TimingDetails | Belos::FinalSummary );
    m_solver_parameter_list->set( "Block Size", 1 );
    m_solver_parameter_list->set( "Num Blocks", 400 );
    m_solver_parameter_list->set( "Maximum Iterations", 500 );
    m_solver_parameter_list->set( "Convergence Tolerance", 1.0e-4 );
    m_solver_parameter_list->set( "Num Recycled Blocks", 300 );

    update_parameters();
  }

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

//.........这里部分代码省略.........
  comm->barrier();
  tm = rcp(new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: 1.5 - MueLu read XML")));
  ParameterListInterpreter mueLuFactory(xmlFileName, *comm);
  comm->barrier();
  tm = Teuchos::null;

  tm = rcp(new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: 2 - MueLu Setup")));

  RCP<Hierarchy> H = mueLuFactory.CreateHierarchy();

  H->GetLevel(0)->Set("A",           A);
  H->GetLevel(0)->Set("Nullspace",   nullspace);

  mueLuFactory.SetupHierarchy(*H);

  comm->barrier();
  tm = Teuchos::null;

  // =========================================================================
  // System solution (Ax = b)
  // =========================================================================

    //
    // generate exact solution using a direct solver
    //
    RCP<Epetra_Vector> exactLsgVec = rcp(new Epetra_Vector(X->Map()));
    {
      fancyout << "========================================================\nCalculate exact solution." << std::endl;
      tm = rcp(new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: 3 - direct solve")));
      exactLsgVec->PutScalar(0.0);
      exactLsgVec->Update(1.0,*X,1.0);
      Epetra_LinearProblem epetraProblem(epA.get(), exactLsgVec.get(), B.get());

      Amesos amesosFactory;
      RCP<Amesos_BaseSolver> rcp_directSolver = Teuchos::rcp(amesosFactory.Create("Amesos_Klu", epetraProblem));
      rcp_directSolver->SymbolicFactorization();
      rcp_directSolver->NumericFactorization();
      rcp_directSolver->Solve();

      comm->barrier();
      tm = Teuchos::null;
    }

    //
    // Solve Ax = b using AMG as a preconditioner in AztecOO
    //

    RCP<Epetra_Vector> precLsgVec = rcp(new Epetra_Vector(X->Map()));
    {
      fancyout << "========================================================\nUse multigrid hierarchy as preconditioner within CG." << std::endl;
      tm = rcp(new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: 4 - AMG as preconditioner")));

      precLsgVec->PutScalar(0.0);
      precLsgVec->Update(1.0,*X,1.0);
      Epetra_LinearProblem epetraProblem(epA.get(), precLsgVec.get(), B.get());

      AztecOO aztecSolver(epetraProblem);
      aztecSolver.SetAztecOption(AZ_solver, AZ_gmres);

      MueLu::EpetraOperator aztecPrec(H);
      aztecSolver.SetPrecOperator(&aztecPrec);

      int maxIts = 50;

      aztecSolver.Iterate(maxIts, tol);

//.........这里部分代码省略.........
    smooProto = Teuchos::rcp( new TrilinosSmoother(ifpackType, ifpackList) );
    RCP<SmootherFactory> SmooFact;
    if (maxLevels > 1)
      SmooFact = rcp( new SmootherFactory(smooProto) );

    // design multigrid hierarchy
    FactoryManager M;
    M.SetFactory("P", PFact);
    M.SetFactory("R", RFact);
    M.SetFactory("Nullspace", PFact);
    M.SetFactory("Smoother", SmooFact);
    M.SetFactory("CoarseSolver", SmooFact);

    H->Setup(M, 0, maxLevels);

    comm->barrier();
    tm = Teuchos::null;

    // =========================================================================
    // System solution (Ax = b)
    // =========================================================================

    //
    // generate exact solution using a direct solver
    //
    RCP<Epetra_Vector> exactLsgVec = rcp(new Epetra_Vector(X->Map()));
    {
      fancyout << "========================================================\nCalculate exact solution." << std::endl;
      tm = rcp(new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: 3 - direct solve")));
      exactLsgVec->PutScalar(0.0);
      exactLsgVec->Update(1.0,*X,1.0);
      Epetra_LinearProblem epetraProblem(epA.get(), exactLsgVec.get(), B.get());

      Amesos amesosFactory;
      RCP<Amesos_BaseSolver> rcp_directSolver = Teuchos::rcp(amesosFactory.Create("Amesos_Klu", epetraProblem));
      rcp_directSolver->SymbolicFactorization();
      rcp_directSolver->NumericFactorization();
      rcp_directSolver->Solve();

      comm->barrier();
      tm = Teuchos::null;
    }

    //
    // Solve Ax = b using AMG as a preconditioner in AztecOO
    //

    RCP<Epetra_Vector> precLsgVec = rcp(new Epetra_Vector(X->Map()));
    {
      fancyout << "========================================================\nUse multigrid hierarchy as preconditioner within CG." << std::endl;
      tm = rcp(new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: 4 - AMG as preconditioner")));

      precLsgVec->PutScalar(0.0);
      precLsgVec->Update(1.0,*X,1.0);
      Epetra_LinearProblem epetraProblem(epA.get(), precLsgVec.get(), B.get());

      AztecOO aztecSolver(epetraProblem);
      aztecSolver.SetAztecOption(AZ_solver, AZ_gmres);

      MueLu::EpetraOperator aztecPrec(H);
      aztecSolver.SetPrecOperator(&aztecPrec);

      int maxIts = 100;
      //double tol2 = 1e-8;

      aztecSolver.Iterate(maxIts, tol);

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

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