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

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

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

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

  // 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

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