pptx - The Trilinos Project

Report
A Tutorial on
Anasazi and Belos
2011 Trilinos User Group Meeting
November 1st, 2011
Chris Baker
David Day
Mike Heroux
Mark Hoemmen
Rich Lehoucq
Mike Parks
Heidi 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.
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
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
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
ARnoldi PACKage
(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
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, matrixvector products, and linear solves
 Allows the creation / registration of new matrix-vector products
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.
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
Belos Solver 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, orthogonalization techniques
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
Anasazi Categories & Solvers
 Anasazi provides solvers for:
 Standard eigenvalue problems: AX = X
 Generalized eigenvalue problems: AX = BX
 Hermitian Eigenproblems
 Block Davidson
 Locally-Optimal Block Preconditioned Conjugate Gradient (LOBPCG)
 Implicit Riemannian Trust-Region (IRTR)
 Non-Hermitian Eigenproblems
 Block Krylov-Schur (BKS)
Anasazi and Belos
(Framework Overview)
GMRES Example
Anasazi and Belos
(Framework Overview)
Multivector
Classes
GMRES Example
Anasazi and Belos
(Framework Overview)
Problem Classes / Operator Classes
Multivector
Classes
GMRES Example
Anasazi and Belos
(Framework Overview)
Problem Classes / Operator Classes
Multivector
Classes
[Mat]OrthoMan
ager
Class
(ICGS, IMGS,
DGKS, SVQB,
TSQR)
GMRES Example
Anasazi and Belos
(Framework Overview)
Problem Classes / Operator Classes
Multivector
Classes
[Mat]OrthoMan
ager
Class
(ICGS, IMGS,
DGKS, SVQB,
TSQR)
StatusTest
Class
GMRES Example
Anasazi and Belos
(Framework Overview)
Problem Classes / Operator Classes
Multivector
Classes
Iteration
Class
[Mat]OrthoMan
ager
Class
(ICGS, IMGS,
DGKS, SVQB,
TSQR)
StatusTest
Class
GMRES Example
Anasazi and Belos
SolverManager Class
(Framework Overview)
Problem Classes / Operator Classes
Multivector
Classes
Iteration
Class
[Mat]OrthoMan
ager
Class
(ICGS, IMGS,
DGKS, SVQB,
TSQR)
StatusTest
Class
GMRES Example
Anasazi and Belos
SolverManager Class
(Framework Overview)
Problem Classes / Operator Classes
Multivector
Classes
Iteration
Class
[Mat]OrthoMan
ager
Class
(ICGS, IMGS,
DGKS, SVQB,
TSQR)
StatusTest
Class
OutputManager
Class
SortManager
Class
Available Solver Components
SNAUPD, DNAUPD, CNAUPD, and ZNAUPD for solving non-Hermit ian e
Separat e subrout ines are required for each of t he four FORT RAN 77
types (single and double precision real, and single and double precis
Moreover, four addit ional subrout ines are needed for a dist ribut ed m
ment at ion. By t emplat ing abst ract numerical int erfaces on operat or
and scalar types, it is only necessary t o maint ain a single code usin
framework.
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
Fig. 1. A n eigensolver
t emplat ed on scalar (ST ), mult ivect or (M V ), and opera
application of an operator
to a multivector.
 Allows user to leverage existing linear algebra software investment.
Anot her aspect of software reuse t hat t emplat ing alleviat es is t hr
 Implementations
arat ion of t he eigensolver algorit hm from t he linear algebra dat a st r
• MultiVecTraits<double,Epetra_MultiVector>
separat ion, as shown in Figure 3.1, allows a user of t he Anasazi frame
age any exist ing linear algebra software invest
• MultiVecTraits<ST,Thyra::MultiVectorBase<ST>
> ment . All t hat is r
t emplat e inst ant iat ion of t he t rait classes, Mul t i VecTr ai>t s and Ope
• MultiVecTraits<ST,Tpetra::MultiVector<ST,LO,GO,Node>
for t he user-defined mult ivect or and operat or, respect ively. T he S
class and respect ive t emplat e inst ant iat ions for different scalar type
Anasazi Eigenproblem Interface
 Provides an interface between the basic iterations and the
eigenproblem to be solved.
 Abstract base class Anasazi::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 class Anasazi::BasicEigenproblem
 Describes standard or general, Hermitian or non-Hermitian
eigenproblems.
Belos LinearProblem Interface
 Provides an interface between the basic iterations and the linear
problem to be solved.
 Templated class Belos::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(…)
Orthogonalization Manager
 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;
int projectAndNormalize(…) const;
magnitudeType orthogError(…) const;
magnitudeType orthonormError(…) const;
 Concrete classes:
Anasazi::
BasicOrthoManager
ICGSOrthoManager
SVQBOrthoManager
Belos::
DGKSOrthoManager
ICGSOrthoManager
IMGSOrthoManager
TSQROrthoManager*
StatusTest Interface




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





StatusType checkStatus( Iteration<…>* solver );
StatusType getStatus() const;
void clearStatus();
void reset();
ostream& print( ostream& os, int indent = 0 ) const;
 StatusType: Passed, Failed, Undefined
 Iteration proceeds until StatusTest returns Passed
 Concrete classes:
Anasazi::
StatusTestMaxIters
StatusTestResNorm
StatusTestOrderedResNorm
StatusTestOutput
Belos::
StatusTestMaxIters
StatusTest[Imp/Gen]ResNorm
StatusTestResNormOutput
StatusTestGeneralOutput
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( int vbLevel );
• int getVerbosity( );
• ostream& stream( MsgType type );
 Query Methods:
• bool isVerbosity( 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.
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
Anasazi Eigensolver Interface
 Provides an abstract interface to Anasazi basic iterations.
 Abstract base class Anasazi::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
Anasazi Eigensolver Manager
 Provides an abstract interface to Anasazi solver managers
(solver strategies)
 Abstract base class Anasazi::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
Belos Iteration Interface


Provides an abstract interface to Belos basic iteration kernels.
Abstract base class Belos::Iteration<ST,MV,OP>







int getNumIters() const; void resetNumIters(int iter);
Teuchos::RCP<const MV> getNativeResiduals( … ) const;
Teuchos::RCP<const MV> getCurrentUpdate() const;
int getBlockSize() const; void setBlockSize(int blockSize);
const LinearProblem<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
Belos Solver Manager

Provides an abstract interface to Belos solver managers (solver strategies)

Abstract base class Belos::SolverManager







void setProblem(…); void setParameters(…);
const Belos::LinearProblem<ST,MV,OP>& getProblem() const;
Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const;
Belos::ReturnType solve();
int getNumIters();
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
Simple Examples
Anasazi Eigensolver Example
(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!!!
…
}
Anasazi Eigensolver Example
(Construct and run the eigensolver)
// Create parameter list to pass into the solver manager
Teuchos::ParameterList MyPL;
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::ReturnType returnCode = MySolverMan.solve();
if (returnCode != Anasazi::Converged) {
// Solver failed!!!
// But wait, we may still have some eigenpairs.
}
Anasazi Eigensolver Example
(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
…
}
…
Belos Linear Solver Example
(Construct the linear problem)
int main(int argc, char *argv[]) {
MPI_Init(&argc,&argv);
Epetra_MpiComm Comm(MPI_COMM_WORLD);
int MyPID = Comm.MyPID();
typedef
typedef
typedef
typedef
typedef
typedef
typedef
double
Teuchos::ScalarTraits<ST>
SCT::magnitudeType
Epetra_MultiVector
Epetra_Operator
Belos::MultiVecTraits<ST,MV>
Belos::OperatorTraits<ST,MV,OP>
ST;
SCT;
MT;
MV;
OP;
MVT;
OPT;
Parameters for
Templates
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);
Get linear
system from
disk
Trilinos/packages/belos/epetra/example/BlockGmres/BlockGmresEpetraExFile.cpp
Belos Linear Solver Example
(Set the solver parameters)
bool verbose = false, debug = false, proc_verbose = false;
int frequency = -1;
// frequency of status test output.
int blocksize = 1;
// blocksize
int numrhs = 1;
// number of right-hand sides to solve for
int maxiters = 100;
// maximum number of iterations allowed
int maxsubspace = 50;
// maximum number of blocks
int maxrestarts = 15;
// number of restarts allowed
MT tol = 1.0e-5;
// relative residual tolerance
Solver
Parameters
const int NumGlobalElements = B->GlobalLength();
ParameterList belosList;
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 );
ParameterList for
SolverManager
Trilinos/packages/belos/epetra/example/BlockGmres/BlockGmresEpetraExFile.cpp
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;
}
LinearProblem
Object
Template Parameters
// 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
SolverManager Object
Belos::ReturnType ret = newSolver->solve();
if (ret!=Belos::Converged) {
std::cout << std::endl << "ERROR:
return -1;
}
std::cout << std::endl << "SUCCESS:
return 0;
Belos did not converge!" << std::endl;
Belos converged!" << std::endl;
Trilinos/packages/belos/epetra/example/BlockGmres/BlockGmresEpetraExFile.cpp
Other Interfaces to
Anasazi / Belos
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]
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");
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 Trilinos Tutorial
http://trilinos.sandia.gov/Trilinos10.8Tutorial.pdf
• See website for more:
http://trilinos.sandia.gov/packages/belos
http://trilinos.sandia.gov/packages/anasazi

similar documents