1 / 42

Anasazi and Belos: Tutorial on Next-Generation Linear Solvers

Explore the capabilities of Anasazi and Belos, the next-generation linear solver libraries, for solving sparse, matrix-free systems. Learn how to use these libraries through simple examples.

kstrine
Download Presentation

Anasazi and Belos: Tutorial on Next-Generation Linear Solvers

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. A Tutorial on Anasazi and Belos2011 Trilinos User Group Meeting November 1st, 2011Chris BakerDavid DayMike HerouxMark HoemmenRich LehoucqMike ParksHeidi Thornquist (Lead) 2011-8264P Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-AC04-94AL85000.

  2. Outline • Belos and Anasazi Framework • Background / Motivation • Framework overview • Available solver components • Using Anasazi and Belos • Simple examples • Through Stratimikos (Belos) • Through LOCA (Anasazi) • Summary

  3. Background / Motivation • Several iterative linear solver / eigensolver libraries exist: • PETSc, SLAP, LINSOL, Aztec(OO), … • SLEPc, PRIMME, ARPACK, … • None of the linear solver libraries can efficiently deal with multiple right-hand sides or sequences of linear systems • Stopping criteria are predetermined for most libraries • The underlying linear algebra is static

  4. AztecOO • A C++ wrapper around Aztec library written in C • Algorithms: GMRES, CG, CGS, BiCGSTAB, TFQMR • Offers status testing capabilities • Output verbosity level can be determined by user • Interface requires Epetra objects • Double-precision arithmetic • Interface to matrix-vector product is defined by the user through the Epetra_Operator

  5. ARnoldiPACKage(ARPACK) Written in Fortran 77 Algorithms: Implicitly Restarted Arnoldi/Lanczos Static convergence tests Output formatting, verbosity level is determined by user Uses LAPACK/BLAS to perform underlying vector space operations Offers abstract interface to matrix-vector products through reverse communication

  6. Scalable Library for Eigenvalue Problem Computations (SLEPc) • Written in C (Campos, Román, Romero, and Thomás, 2011). • Provides some basic eigensolvers as well as wrappers around: • ARPACK(Lehoucq, Maschhoff, Sorensen, and Yang, 1998) • PRIMME (Stathopoulos, 2006) • BLOPEX (Knyazev, 2007) • BLZPACK (Marques, 1995) • TRLAN(Wu and Simon, 2001) • Native Algorithms: Power/Subspace Iteration, RQI, Arnoldi, KS • Wrapped Algorithms: IRAM/IRLM (ARPACK), Block Lanczos (BLZPACK), LOBPCG (BLOPEX), and Lanczos(TRLAN) • Static convergence tests • Uses PETSc to perform underlying vector space operations, matrix-vector products, and linear solves • Allows the creation / registration of new matrix-vector products

  7. Anasazi and Belos • Next generation linear solver (Belos) and eigensolver (Anasazi) libraries, written in templated C++. • Iterative methods for solving sparse, matrix-free systems • Provide a generic interface to a collection of algorithms for solving linear problems and eigenproblems. • Algorithms developed with generic programming techniques. • Algorithmic components: • Ease the implementation of complex algorithms • Operator/MultiVector interface (and Teuchos::ScalarTraits): • Allow the user to leverage their existing software investment • Multi-precision solver capability • Design offers: Interoperability, extensibility, and reusability • Includes block linear solvers and eigensolvers.

  8. Why are Block Solvers Useful? • In general, block solvers enable the use of faster computational kernels. • Block Eigensolvers ( Op(A)X = X ): • Reliably determine multiple and/or clustered eigenvalues. • Example applications: • Stability analysis / Modal analysis • Bifurcation analysis (LOCA) • Block Linear Solvers ( Op(A)X = B ): • Useful for when multiple solutions are required for the same system of equations. • Example applications: • Perturbation analysis • Optimization problems • Single right-hand sides where A has a handful of small eigenvalues • Inner-iteration of block eigensolvers

  9. BelosSolver Categories • Belos provides solvers for: • Single RHS: Ax = b • Multiple RHS (available simultaneously): AX = B • Multiple RHS (available sequentially): Axi = bi , i=1,…,k • Sequential Linear systems: Aixi = bi , i=1,…,k • Linear Least Squares: min || Ax – b ||2 • Leverage research advances of solver community: • Block methods: block GMRES [Vital], block CG/BICG [O’Leary] • “Seed” solvers: hybrid GMRES [Nachtigal, et al.] • “Recycling” solvers for sequences of linear systems [Parks, et al.] • Restarting, orthogonalizationtechniques

  10. Belos Solvers • Hermitian Systems (A = AH) • CG / Block CG • Pseudo-Block CG (Perform single-vector algorithm simultaneously) • RCG (Recycling Conjugate Gradients) • PCPG (Projected CG) • MINRES • Non-Hermitian System (A ≠ AH) • Block GMRES • Pseudo-Block GMRES (Perform single-vector algorithm simultaneously) • Block FGMRES (Variable preconditioner) • Hybrid GMRES • TFQMR • GCRODR / Block GCRODR (Recycling GMRES) • Linear Least Squares • LSQR

  11. AnasaziCategories & Solvers • Anasazi provides solvers for: • Standard eigenvalue problems: AX= X • Generalized eigenvalue problems: AX = BX • HermitianEigenproblems • Block Davidson • Locally-Optimal Block Preconditioned Conjugate Gradient (LOBPCG) • Implicit Riemannian Trust-Region (IRTR) • Non-HermitianEigenproblems • Block Krylov-Schur (BKS)

  12. Anasazi and Belos(Framework Overview) GMRES Example

  13. Anasazi and Belos(Framework Overview) Multivector Classes GMRES Example

  14. Anasazi and Belos(Framework Overview) Problem Classes / Operator Classes Multivector Classes GMRES Example

  15. Anasazi and Belos(Framework Overview) Problem Classes / Operator Classes Multivector Classes [Mat]OrthoManagerClass (ICGS, IMGS, DGKS, SVQB, TSQR) GMRES Example

  16. Anasazi and Belos(Framework Overview) Problem Classes / Operator Classes Multivector Classes [Mat]OrthoManagerClass (ICGS, IMGS, DGKS, SVQB, TSQR) StatusTest Class GMRES Example

  17. Anasazi and Belos(Framework Overview) Problem Classes / Operator Classes Iteration Class Multivector Classes [Mat]OrthoManagerClass (ICGS, IMGS, DGKS, SVQB, TSQR) StatusTest Class GMRES Example

  18. Anasazi and Belos(Framework Overview) SolverManager Class Problem Classes / Operator Classes Iteration Class Multivector Classes [Mat]OrthoManagerClass (ICGS, IMGS, DGKS, SVQB, TSQR) StatusTest Class GMRES Example

  19. Anasazi and Belos(Framework Overview) SolverManager Class Problem Classes / Operator Classes Iteration Class Multivector Classes [Mat]OrthoManagerClass (ICGS, IMGS, DGKS, SVQB, TSQR) StatusTest Class OutputManager Class SortManager Class

  20. Available Solver Components

  21. Linear Algebra Interface • MultiVecTraits • Abstract interface to define the linear algebra required by most iterative solvers: • creational methods • dot products, norms • update methods • initialize / randomize • OperatorTraits • Abstract interface to enable the application of an operator to a multivector. • Allows user to leverage existing linear algebra software investment. • Implementations • MultiVecTraits<double,Epetra_MultiVector> • MultiVecTraits<ST,Thyra::MultiVectorBase<ST> > • MultiVecTraits<ST,Tpetra::MultiVector<ST,LO,GO,Node> >

  22. Anasazi Eigenproblem Interface • Provides an interface between the basic iterations and the eigenproblem to be solved. • Abstract base classAnasazi::Eigenproblem • Allows spectral transformations to be removed from the algorithm. • Differentiates between standard and generalized eigenproblems. • Stores number of requested eigenvalues, symmetry, initial vector, auxiliary vectors, stiffness/mass matrix, operator, and eigensolution. • Evecs, Evals, Espace, index, numVecs • Concrete classAnasazi::BasicEigenproblem • Describes standard or general, Hermitian or non-Hermitianeigenproblems.

  23. BelosLinearProblemInterface • Provides an interface between the basic iterations and the linear problem to be solved. • Templated classBelos::LinearProblem<ST,MV,OP> • Allows preconditioning to be removed from the algorithm. • Behavior defined through traits mechanisms. • Methods: • setOperator(…) / getOperator() • setLHS(…) / getLHS() • setRHS(…) / getRHS() • setLeftPrec(…) / getLeftPrec() / isLeftPrec() • setRightPrec(…) / getRightPrec() / isRightPrec() • apply(…) / applyOp(…) / applyLeftPrec(…) / applyRightPrec(…) • setHermitian(…) / isHermitian() • setProblem(…)

  24. OrthogonalizationManager • Abstract interface to orthogonalization / orthonormalization routines for multivectors. • Abstract base class [Anasazi/Belos]::[Mat]OrthoManager • void innerProd(…) const; • void norm(…) const; • void project(…) const; • int normalize(…) const; • intprojectAndNormalize(…) const; • magnitudeTypeorthogError(…) const; • magnitudeTypeorthonormError(…) const; • Concrete classes: Anasazi:: BasicOrthoManager ICGSOrthoManager SVQBOrthoManager Belos:: DGKSOrthoManager ICGSOrthoManager IMGSOrthoManager TSQROrthoManager*

  25. StatusTestInterface • Informs solver iterate of change in state, as defined by user. • Similar to NOX / AztecOO. • Multiple criterion can be logically connected using StatusTestCombo • Abstract base class [Anasazi/Belos]::StatusTest • StatusTypecheckStatus( Iteration<…>* solver ); • StatusTypegetStatus() const; • void clearStatus(); • void reset(); • ostream& print( ostream& os, int indent = 0 ) const; • StatusType: Passed, Failed, Undefined • Iteration proceeds untilStatusTestreturnsPassed • Concrete classes: Anasazi:: StatusTestMaxIters StatusTestResNorm StatusTestOrderedResNorm StatusTestOutput Belos:: StatusTestMaxIters StatusTest[Imp/Gen]ResNorm StatusTestResNormOutput StatusTestGeneralOutput

  26. Output Manager Interface • Templated class that enables control of the linear solver output. • Behavior defined through traits mechanism • OutputManager<ST> • Get/Set Methods: • void setVerbosity( intvbLevel ); • intgetVerbosity( ); • ostream& stream( MsgType type ); • Query Methods: • boolisVerbosity( MsgType type ); • Print Methods: • void print( MsgType type, const string output ); • Message Types: • Errors, Warnings, IterationDetails, OrthoDetails, FinalSummary, TimingDetails, StatusTestDetails, Debug • Default is lowest verbosity (Errors), output on one processor.

  27. Anasazi Sort Manager • Abstract interface for managing the sorting of the eigenvalues computed by the eigensolver. • Important tool to complement spectral transformations. • Only two methods: • ReturnType sort(Eigensolver<ST,MV,OP>* solver, int n, ST *evals, std::vector<int> *perm=0); • ReturnType sort(Eigensolver<ST,MV,OP>* solver, int n, ST *r_evals, ST *i_evals,std::vector<int> *perm=0); • Concrete class Anasazi::BasicSort • Provides basic sorting methods: • largest/smallest magnitude • largest/smallest real part • largest/smallest imaginary part

  28. Anasazi Eigensolver Interface • Provides an abstract interface to Anasazi basic iterations. • Abstract base classAnasazi::Eigensolver • get / reset number of iterations • native residuals • current / maximum subspace size • set / get auxiliary vectors • eigenproblem • initialize / iterate • Concrete solvers (iterations): • Anasazi::BlockKrylovSchur • Anasazi::BlockDavidson • Anasazi::LOBPCG • Anasazi::IRTR

  29. Anasazi Eigensolver Manager • Provides an abstract interface to Anasazi solver managers (solver strategies) • Abstract base classAnasazi::SolverManager • Access to the eigenproblem • Solve the eigenproblem • Solvers are parameter list driven, take two arguments: • Anasazi::Eigenproblem • Teuchos::ParameterList • Concrete solver managers: • Anasazi::BlockKrylovSchurSolMgr • Anasazi::BlockDavidsonSolMgr • Anasazi::LOBPCGSolMgr • Anasazi::IRTRSolMgr

  30. Belos Iteration Interface • Provides an abstract interface to Belos basic iteration kernels. • Abstract base class Belos::Iteration<ST,MV,OP> • intgetNumIters() const; void resetNumIters(intiter); • Teuchos::RCP<const MV> getNativeResiduals( … ) const; • Teuchos::RCP<const MV> getCurrentUpdate() const; • intgetBlockSize() const; void setBlockSize(intblockSize); • constLinearProblem<ST,MV,OP>& getProblem() const; • void iterate(); void initialize(); • Iterations require these classes: • Belos::LinearProblem, Belos::OutputManager, Belos::StatusTest, Belos::OrthoManager • Implementations: • Belos::BlockGmresIter • Belos::BlockFGmresIter • Belos::PseudoBlock[CG/Gmres]Iter • Belos::[Block]GCRODRIter • Belos::[CG/BlockCG]Iter, • Belos::RCGIter • Belos::PCPGIter • Belos::TFQMRIter • Belos::LSQRIter • Belos::MinresIter

  31. Belos Solver Manager • Provides an abstract interface to Belos solver managers (solver strategies) • Abstract base class Belos::SolverManager • void setProblem(…); void setParameters(…); • constBelos::LinearProblem<ST,MV,OP>& getProblem() const; • Teuchos::RCP<constTeuchos::ParameterList> getValidParameters() const; • Teuchos::RCP<constTeuchos::ParameterList> getCurrentParameters() const; • Belos::ReturnType solve(); • intgetNumIters(); • Solvers are parameter list driven, take two arguments: • Belos::LinearProblem • Teuchos::ParameterList[validated] • Implementations: • Belos::BlockGmresSolMgr • Belos::PseudoBlock[CG/Gmres]SolMgr • Belos::[Block]GCRODRSolMgr • Belos::BlockCGSolMgr • Belos::RCGSolMgr • Belos::PCPGSolMgr • Belos::TFQMRSolMgr • Belos::LSQRSolMgr • Belos::MinresSolMgr

  32. Simple Examples

  33. Anasazi EigensolverExample(Construct the eigenproblem) // Create eigenproblem const int nev = 5; RefCountPtr<Anasazi::BasicEigenproblem<ScalarType,MV,OP> > problem = Teuchos::rcp( new Anasazi::BasicEigenproblem<ScalarType,MV,OP>(K,M,ivec) ); // // Inform the eigenproblem that it is Hermitian problem->setHermitian(true); // // Set the number of eigenvalues requested problem->setNEV( nev ); // // Inform the eigenproblem that you are done passing it information bool ret = problem->setProblem(); if (boolret != true) { // Anasazi::BasicEigenproblem::SetProblem() returned with error!!! … }

  34. Anasazi EigensolverExample(Construct and run the eigensolver) // Create parameter list to pass into the solver manager Teuchos::ParameterListMyPL; MyPL.set( "Verbosity", Anasazi::Errors + Anasazi::Warnings ); MyPL.set( "Which", "SM"); MyPL.set( "Block Size", 5 ); MyPL.set( "Num Blocks", 8 ); MyPL.set( "Maximum Restarts", 100 ); MyPL.set( "Convergence Tolerance", 1.0e-6 ); MyPL.set( "Use Locking", true ); MyPL.set( "Locking Tolerance", tol/10 ); // // Create the solver manager Anasazi::BlockDavidsonSolMgr<ScalarType,MV,OP> MySolverMan(problem, MyPL); // Solve the problem to the specified tolerances or length Anasazi::ReturnTypereturnCode = MySolverMan.solve(); if (returnCode != Anasazi::Converged) { // Solver failed!!! // But wait, we may still have some eigenpairs. }

  35. Anasazi EigensolverExample(Retrieve the eigenpairs) // Get the eigenvalues and eigenvectors from the eigenproblem Anasazi::Eigensolution<ScalarType,MV> sol = problem->getSolution(); std::vector<MagnitudeType> evals = sol.Evals; RefCountPtr<MV> evecs = sol.Evecs; int numev = sol.numVecs; // // Check to see if there are any computed eigenpairs if (numev > 0) { // Do something with computed eigenpairs … } …

  36. Belos Linear Solver Example(Construct the linear problem) int main(intargc, char *argv[]) { MPI_Init(&argc,&argv); Epetra_MpiCommComm(MPI_COMM_WORLD); intMyPID = Comm.MyPID(); typedef double ST; typedefTeuchos::ScalarTraits<ST> SCT; typedef SCT::magnitudeType MT; typedefEpetra_MultiVector MV; typedefEpetra_Operator OP; typedefBelos::MultiVecTraits<ST,MV> MVT; typedefBelos::OperatorTraits<ST,MV,OP> OPT; using Teuchos::ParameterList; using Teuchos::RCP; using Teuchos::rcp; // Get the problem std::string filename("orsirr1.hb"); RCP<Epetra_Map> Map; RCP<Epetra_CrsMatrix> A; RCP<Epetra_MultiVector> B, X; RCP<Epetra_Vector> vecB, vecX; EpetraExt::readEpetraLinearSystem(filename, Comm, &A, &Map, &vecX, &vecB); X = Teuchos::rcp_implicit_cast<Epetra_MultiVector>(vecX); B = Teuchos::rcp_implicit_cast<Epetra_MultiVector>(vecB); Parameters for Templates Get linear system from disk Trilinos/packages/belos/epetra/example/BlockGmres/BlockGmresEpetraExFile.cpp

  37. Belos Linear Solver Example(Set the solver parameters) bool verbose = false, debug = false, proc_verbose = false; int frequency = -1; // frequency of status test output. intblocksize = 1; // blocksize intnumrhs = 1; // number of right-hand sides to solve for intmaxiters = 100; // maximum number of iterations allowed intmaxsubspace = 50; // maximum number of blocks intmaxrestarts = 15; // number of restarts allowed MT tol = 1.0e-5; // relative residual tolerance const intNumGlobalElements = B->GlobalLength(); ParameterListbelosList; belosList.set( "Num Blocks", maxsubspace); // Maximum number of blocks in Krylov factorization belosList.set( "Block Size", blocksize ); // Blocksize to be used by iterative solver belosList.set( "Maximum Iterations", maxiters ); // Maximum number of iterations allowed belosList.set( "Maximum Restarts", maxrestarts ); // Maximum number of restarts allowed belosList.set( "Convergence Tolerance", tol ); // Relative convergence tolerance requested int verbosity = Belos::Errors + Belos::Warnings; if (verbose) { verbosity += Belos::TimingDetails + Belos::StatusTestDetails; if (frequency > 0) belosList.set( "Output Frequency", frequency ); } if (debug) { verbosity += Belos::Debug; } belosList.set( "Verbosity", verbosity ); Solver Parameters ParameterList for SolverManager Trilinos/packages/belos/epetra/example/BlockGmres/BlockGmresEpetraExFile.cpp

  38. Belos Linear Solver Example(Solve the linear problem) // Construct linear problem instance. Belos::LinearProblem<double,MV,OP> problem( A, X, B ); bool set = problem.setProblem(); if (set == false) { std::cout << std::endl << "ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl; return -1; } // Start block GMRES iteration Belos::OutputManager<double> My_OM(); // Create solver manager. RCP< Belos::SolverManager<double,MV,OP> > newSolver = rcp( new Belos::BlockGmresSolMgr<double,MV,OP>(rcp(&problem,false), rcp(&belosList,false))); // Solve Belos::ReturnType ret = newSolver->solve(); if (ret!=Belos::Converged) { std::cout << std::endl << "ERROR: Belos did not converge!" << std::endl; return -1; } std::cout << std::endl << "SUCCESS: Belos converged!" << std::endl; return 0; LinearProblem Object Template Parameters SolverManager Object Trilinos/packages/belos/epetra/example/BlockGmres/BlockGmresEpetraExFile.cpp

  39. Other Interfaces to Anasazi / Belos

  40. Stratimikos-Belos Interface • Stratimikos: • Thyra-based linear solver and preconditioner strategy package • MultiVecTraits<ST,Thyra::MultiVectorBase<ST> > • OperatorTraits<ST,Thyra::MultiVectorBase<ST>,Thyra::LinearOpBase<ST> > • Stratimikos-Belos Interface • Intent: access current Belos solver managers using a factory interface • Implements Thyra::LinearOpWithSolveFactory / Thyra::LinearOpWithSolve • Stratimikos interface uses valid parameter list generated from Belos • Check out: stratimikos/adapters/belos/example/LOWSFactory/[Epetra/Tpetra]

  41. Linear Stability Analysis Through Anasazi and LOCA • LOCA provides parameter continuation and bifurcation tracking for large-scale codes • Changes in stability of steady-states indicated by eigenvalues of linearized system crossing imaginary axis • LOCA provides interface to Anasazi to compute eigenvalues • Interfaced through NOX/LOCA abstract layers • No additional work necessary once LOCA is supported • LOCA provides spectral transformations to emphasize eigenvalues near imaginary axis • Jacobian-inverse – • Shift-invert – • Cayley – stepperList.set("Compute Eigenvalues", true); Teuchos::ParameterList& aList = stepperList.sublist("Eigensolver"); aList.set("Method", "Anasazi"); aList.set("Block Size", 1) aList.set("Num Blocks", 50); aList.set("Num Eigenvalues", 3) aList.set("Step Size", 1); aList.set("Maximum Restarts",1); aList.set("Operator", "Cayley"); aList.set("Cayley Pole", 0.1); aList.set("Cayley Zero", -0.1); aList.set("Sorting Order", "CA");

  42. Summary • Belos and Anasazi are next-generation linear and eigensolver libraries • Designed for interoperability, extensibility, and reusability • Belos and Anasazi are readily available: • Can be used as standalone linear and eigensolvers • Belos available through Stratimikos • Anasazi available through LOCA • Check out the TrilinosTutorial http://trilinos.sandia.gov/Trilinos10.8Tutorial.pdf • See website for more: http://trilinos.sandia.gov/packages/belos http://trilinos.sandia.gov/packages/anasazi

More Related