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) 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.

2011-8264P

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, matrix vector 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 Sol ver Categories

 Belos provides solvers for:      Single RHS:

Ax = b

Multiple RHS (available simultaneously):

AX = B

Multiple RHS (available sequentially):

Ax i

Sequential Linear systems:

A i x i = b i , i =1,…,k = b i , 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 = A H )      CG / Block CG Pseudo-Block CG (Perform single-vector algorithm simultaneously) RCG (Recycling Conjugate Gradients) PCPG (Projected CG) MINRES Non-Hermitian System (A ≠ A H )       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

Multivector Classes

Anasazi and Belos

(Framework Overview) GMRES Example

Multivector Classes

Anasazi and Belos

(Framework Overview) Problem Classes / Operator Classes GMRES Example

Multivector Classes

Anasazi and Belos

(Framework Overview) Problem Classes / Operator Classes [Mat]OrthoMan ager Class (ICGS, IMGS, DGKS, SVQB, TSQR) GMRES Example

Multivector Classes StatusTest Class

Anasazi and Belos

(Framework Overview) Problem Classes / Operator Classes [Mat]OrthoMan ager Class (ICGS, IMGS, DGKS, SVQB, TSQR) GMRES Example

Multivector Classes StatusTest Class

Anasazi and Belos

(Framework Overview) Problem Classes / Operator Classes Iteration Class [Mat]OrthoMan ager Class (ICGS, IMGS, DGKS, SVQB, TSQR) GMRES Example

Anasazi and Belos

(Framework Overview) SolverManager Class Problem Classes / Operator Classes Iteration Class Multivector Classes StatusTest Class [Mat]OrthoMan ager Class (ICGS, IMGS, DGKS, SVQB, TSQR) GMRES Example

Anasazi and Belos

(Framework Overview) SolverManager Class Problem Classes / Operator Classes Iteration Class Multivector Classes StatusTest Class [Mat]OrthoMan ager Class (ICGS, IMGS, DGKS, SVQB, TSQR) OutputManager Class SortManager Class

Available Solver Components

8

·

C. G. Baker et al   We not e t hat t he Anasazi framework was designed t o support block met hods, defined as those that apply A or B to a collection of vectors (a multivector). One advant age of using a mult ivect or dat a st ruct ure is t o improve t he rat io of float ing point operat ions t o memory references and so bet t er exploit a memory hierarchy.

Templat ing an eigensolver on operat or, mult ivect or, and scalar types makes soft ware reuse easier. Consider in cont rast t hat ARPACK implement s t he subrout ines SNAUPD, DNAUPD, CNAUPD, and ZNAUPD for solving non-Hermit ian eigenproblems.

Separat e subrout ines are required for each of t he four FORT RAN 77 floating point types (single and double precision real, and single and double precision complex).

Linear Algebra Interface

ment at ion. By t emplat ing abst ract numerical int erfaces on operat or, mult ivect or,

MultiVecTraits

and scalar types, it is only necessary t o maint ain a single code using t he Anasazi framework.

 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   Allows user to leverage existing linear algebra software investment.

Anot her aspect of software reuse t hat t emplat i ng alleviat es is t hrough t he sep Implementations arat ion of t he eigensolver algorit hm from t he linear algebra dat a st ruct ures. T his • • •

MultiVecTraits >

for t he us er-defined multivector and operator, respectively. The Scal ar Tr ai t s class and respect ive t emplat e inst ant iat ions for different scalar types are provided by t he Trilinos Teuchos package [Heroux et al. ]. Anot her friendly aspect of employ ing t emplat es and t rait s mechanisms is t hat t he Anasazi eigensolver, eigenpr oblem, and eigensolut ion are all defined by the specified scalar, multivect or, and opera t or type at compile t ime. T his approach, as opposed t o using abst ract int erfaces and dynamic polymorphism, avoids any dynamic cast ing of t he mult ivect ors and operat ors in t he us er’s interaction with the Anasazi framework.

T he Mul t i VecTr ai t s and Oper at or Tr ai t s classes specify t he operat ions t hat t he mult ivect or and operat or type must support in order for t hem t o be used by Anasazi.

T hrough t he observat ions made in Sect ion 2, it is clear t hat t he Oper at or Tr ai t s class only needs t o provide one met hod, described in Table I, t hat applies an operat or t o a mult ivect or. T his int erface defines the only interaction required from an operat or, even t hough t he underlying operat or may be a mat rix, spect ral t ransformat ion, or precondit ioner.

A CM T r ansact i ons on M at hem at i cal Sof t war e, Vol . V , N o. N , Januar y 2008.

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

   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  Concrete classes:

Anasazi:: Belos:: Passed StatusTestMaxIters StatusTestResNorm StatusTestOrderedResNorm StatusTestOutput 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

 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* solver, int n, ST *evals, std::vector *perm=0);



ReturnType sort(Eigensolver* solver, int n, ST *r_evals, ST *i_evals, std::vector *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       get / reset number of iterations native residuals

Anasazi::Eigensolver

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



int getNumIters() const; void resetNumIters(int iter);

    

Teuchos::RCP getNativeResiduals( … ) const; Teuchos::RCP getCurrentUpdate() const; int getBlockSize() const; void setBlockSize(int blockSize); const LinearProblem& 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& getProblem() const; Teuchos::RCP getValidParameters() const;

  

Teuchos::RCP 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 > problem = Teuchos::rcp( new Anasazi::BasicEigenproblem(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 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 sol = problem->getSolution(); std::vector evals = sol.Evals; RefCountPtr 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 double ST; typedef Teuchos::ScalarTraits SCT; typedef SCT::magnitudeType MT; typedef Epetra_MultiVector typedef Epetra_Operator MV; OP; typedef Belos::MultiVecTraits MVT; typedef Belos::OperatorTraits OPT;

Parameters for Templates

using Teuchos::ParameterList; using Teuchos::RCP; using Teuchos::rcp; // Get the problem std::string filename("orsirr1.hb"); RCP Map; RCP A; RCP B, X; RCP vecB, vecX; EpetraExt::readEpetraLinearSystem(filename, Comm, &A, &Map, &vecX, &vecB); X = Teuchos::rcp_implicit_cast(vecX); B = Teuchos::rcp_implicit_cast(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 problem( A, X, B ); bool set = problem.setProblem();

LinearProblem

if (set == false) { set up correctly!" << std::endl;

Object

return -1;

Template Parameters

// Start block GMRES iteration Belos::OutputManager My_OM(); // Create solver manager.

RCP< Belos::SolverManager > newSolver =

SolverManager Object

rcp( new Belos::BlockGmresSolMgr(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; 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 >

•

OperatorTraits,Thyra::LinearOpBase >

 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