# Presentation: 001002

## An Overview of AnasaziNovember 2nd, 2:15-3:15pmChris Baker Ulrich HetmaniukRich LehoucqHeidi Thornquist (Lead)

- Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company,for the United States Department of Energy under contract DE-AC04-94AL85000.

**Notes:**
Questions from audience

Can you use another matrix-vector multiply?
Are multivectors supported?
Are Petsc and Hypre available via TSF?
How does Trilinos differ from Petsc?

## Outline

- Background
- ARPACK
- SLEPc
- Eigensolver design goals
- Generic programming
- Anasazi
- What’s in Trilinos 6.0
- Current Interfaces
- Current Customers
- Future Development
- Concluding remarks

## Background

- Several iterative eigensolver libraries exist:
- SLEPc
- LOBPCG (hypre)
- ARPACK
- ...
- None are (completely) written in C++.
- Stopping criteria are predetermined for most libraries.
- The underlying linear algebra is static.

## ARnoldi PACKage (ARPACK)

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

**Scalable Library for Eigenvalue Problem Computations****  ****( SLEPc )**

- Written in C (Hernández, Román, and Vidal, 2003).
- Provides some basic eigensolvers as well as wrappers around:
- ARPACK (Lehoucq, Maschhoff, Sorensen, and Yang, 1998)
- BLZPACK (Marques, 1995)
- PLANSO (Wu and Simon 1997)
- TRLAN (Wu and Simon, 2001)
- Native Algorithms:  Power/Subspace Iteration, RQI, Arnoldi
- Wrapped Algorithms:  IRAM/IRLM (ARPACK), Block Lanczos   (BLZPACK), and Lanczos (PLANSO/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.

## Eigensolver Software Design

- Provides an abstract interface to an operator-vector products, scaling, and preconditioning.
- Allows the user to enlist any linear algebra package for the elementary vector space operations essential to the algorithm (Epetra, PETSc, etc.).
- Allows the user to define convergence of any algorithm (a.k.a. status testing).
- Allows the user to determine the verbosity level, formatting, and processor for the output.
- Allows these decisions to be made at runtime.

**Eigensolver Software Design**

## AlgorithmA( ... )
	while ( myStatusTest.CheckStatus(...)
			== Unconverged )
     ...
     ...
     % Compute operator-vector product 
     myProblem.ApplyOp(NOTRANS,v,w); 
      = w.dot(v);
     ...
	end
	myOM.print(something);
  	return (solution);
end

**AlgorithmA( ... )**
**	****while ( ****myStatusTest.CheckStatus(****...****)**
**	****	****	****== Unconverged )**
**     ...**
**     ...**
**     ****% Compute operator-vector product **
**     ****myProblem.****ApplyOp(****NOTRANS,****v,w****)****; **
**     ******** ****= ****w****.dot(****v****)****;**
**     ...**
**	****end**
**	****myOM.print(****something****)****;**
**  ****	****return (****solution****);**
**end**

**Generic Eigensolver**

## Benefits of Generic Programming

**Benefits of Generic Programming**

- Generic programming techniques ease the implementation of complex algorithms.
- Developing algorithms with generic programming techniques is easier on the developer, while still allowing them to build a flexible and powerful software package.
- Generic programming techniques also allow the user to leverage their existing software investment.

_**Caveat:**_**  It’s not as easy as taking a piece of code and adding:**
**template <typename OT, typename ST>**
**More work has to be done to handle “numeric traits”**

## Introducing Anasazi

- Next generation eigensolvers library, written in templated C++.
- Provide a generic interface to a collection of algorithms for solving large-scale eigenproblems.
- Algorithm implementation is accomplished through the use of traits clases and abstract base classes:
- e.g.: MultiVecTraits, OperatorTraits
- e.g.: Eigensolver, Eigenproblem, StatusTest
- Includes block eigensolvers:
- Higher operator performance.
- More reliable.
- Solves *AX = X*** or *AX = BX****.*

## Anasazi Status(Trilinos Release 6.0.5)

- Anasazi
- Currently has three solvers:
- Block Krylov-Schur Method
- Block Davidson
- LOBPCG
- Can solve standard and generalized eigenproblems.
- Can solve Hermitian and non-Hermitian eigenproblems.
- Can target largest or smallest eigenvalues.
- Linear algebra adapters for Epetra, NOX/LOCA and Thyra.
- Epetra interface accepts Epetra_Operators, so it can be used with Belos, AztecOO, etc...
- Block size is independent of number of requested eigenvalues.

## Anasazi Design

- Eigenproblem Class
- Describes the problem and stores the answer.
- Eigensolver Class
- Provide iteration interface.
- MultiVecTraits and OperatorTraits
- Traits classes for interfacing linear algebra.
- SortManager, OrthoManager
- Specify these operations.
- StatusTest Class
- Control testing of convergence, etc.
- OutputManager Class
- Control verbosity and printing in a MP scenario.

## Anasazi Eigenproblem Interface

**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.
- Defines the inner product to be used for the problem.
- Differentiates between standard and generalized eigenproblems.
- Initial vector, stiffness/mass matrix, operator, eigenvalues, eigenvectors.
- Concrete class Anasazi::BasicEigenproblem
- Describes standard or general, Hermitian or non-Hermitian eigenproblems and the Euclidean inner product.

## Anasazi Eigensolver Interface

- Provides an abstract interface to Anasazi basic iterations.
- Abstract base class Anasazi::Eigensolver
- number of iterations / restarts
- native residuals
- eigenproblem
- initialize/finalize
- iterate()
- Specialization for solvers that can provide more information
- Block Krylov-Schur -> Ritz values / residuals

## Anasazi Linear Algebra Interface

**Anasazi Linear Algebra Interface**

- Anasazi::MultiVecTraits
- Abstract interface to define the linear algebra required by most iterative eigensolvers:
- creational methods
- dot products, norms
- update methods
- initialize / randomize

- Anasazi::OperatorTraits
- Abstract interface to enable the application of an operator to a multivector.

## Three approaches for using your own linear algebra:
Specialize the traits classes for the multivector and operator:
Fill out Anasazi::MultiVecTraits<ST,MyMV> 
Fill out Anasazi::OperatorTraits<ST,MyMV,MyOP>
Subclass abstract base classes, for which traits specializations are already defined:
class MyMV : public Anasazi::MultiVec<ST>
class MyOP : public Anasazi::Operator<ST>
Implement your linear algebra in the Thyra framework:
Then use the Anasazi adapter to Thyra.

- Three approaches for using your own linear algebra:
- Specialize the traits classes for the multivector and operator:
- Fill out Anasazi::MultiVecTraits<ST,MyMV>
- Fill out Anasazi::OperatorTraits<ST,MyMV,MyOP>
- Subclass abstract base classes, for which traits specializations are already defined:
- class MyMV : public Anasazi::MultiVec<ST>
- class MyOP : public Anasazi::Operator<ST>
- Implement your linear algebra in the Thyra framework:
- Then use the Anasazi adapter to Thyra.

**Interfacing Your LA With Anasazi**

- Or you can simply use Epetra!

## Testing Your Interface

- We need a way to test functionality of a new linear algebra library.
- Anasazi includes two routines for testing LA libraries/interfaces:
- Anasazi::TestMultiVecTraits<ST,MV>(...)
- Anasazi::TestOperatorTraits<ST,MV,OP>(...)
- Tests MultiVecTraits/OperatorTraits specialization and underlying classes.
- Interfaces don’t have explicit access to object data, so tests are limited.

## Anasazi StatusTest Interface

- Allows user to decide when solver is finished.
- Similar to NOX / AztecOO.
- Multiple criterion can be logically connected.
- Abstract base class methods:
- StatusType CheckStatus( Eigensolver<...>* solver );
- StatusType GetStatus() const;
- void Reset();
- ostream& Print( ostream& os ) const;
- Concrete classes:
- Maximum iterations/restarts,
- Residual norm,
- Combinations {AND,OR}, ...

## Anasazi Output Manager Interface

**Anasazi Output Manager Interface**

- Allows user to control which processor will handle output from the solver and the verbosity.
- Default is lowest verbosity, outputting on proc 0.
- Methods:
- Get/Set Methods:
- void SetOStream( ostream& os );
- void SetVerbosity( int vbLevel );
- ostream& GetOStream();
- Query Methods:
- bool isVerbosity( MsgType type );
- bool isVerbosityAndPrint( MsgType type );
- bool doPrint( MsgType type );

## Anasazi Sort Manager Interface

**Anasazi Sort Manager Interface**

- 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

## Current Anasazi Customers

- LOCA
- Bifurcation analysis
- Continuation methods
- RBGen
- Generating bases for model reduction

## Infomercial: RBGen

**Infomercial: RBGen**

- RBGen - **R**educed **B**asis **Gen**erator
- Lead developer: Heidi Thornquist
- Implemented tool to compute reduced basis (RBGen)
- Modular software design facilitates research of RBMs and snapshot generation, Exodus I/O utilities (Paul Lin).
- General purpose tool for any physics application.
- Quick development using tools already written (Anasazi/Teuchos)
- First reduced-basis method: POD
- Centroidal Voronoi Tessellations (CVT)

## Reduced Order Modeling

**Problem:**  Computing high-fidelity state solutions are extremely computationally intensive when performing real-time analysis.
**Solution:  **Reduce the cost of the nonlinear state solution by using a reduced-order model for the state.
- Need for reduced order modeling capabilities:
- Forward problem (flow modeling)
- Optimization problem (source inversion)
- Proper Orthogonal Decomposition (POD)
- Attempt to capture the dominant behavior of the method.
- Constructs basis from a set of solution snapshots.
- Compute first *n *left singular vectors of snapshot matrix.
- Work with John Shadid, Paul Lin, Rich Lehoucq, and Max Gunzburger (Florida State University)

## RBGen::POD

- Given matrix *A* of *n* snapshots of length *m, *we need the first *k* left singular vectors of *A*
- These are the rightmost *k* eigenvectors of *A**H**·A*
- Flexibility of Anasazi allows trivial usage in RBGen:
- Anasazi operator consists of two multiplications by snapshot matrix.
- RBGen has access to any Anasazi eigensolver!

## ROM Preliminary Qualitative Results Transient LES-k 3D Airport Simulation

**Velocity field magnitude averaged from 2200-7200 sec**

**Velocity field magnitude averaged from 2200-7200 sec**

**3D airport (230K nodes; 1.15M unknowns)**
**Large output files (5 unknowns: Vx, Vy, Vz, P, TKE)**
- output every time step for 6000 sec (each time step roughly 1 sec) gives 5.5 GB exodusII file
**7200 seconds of simulation time required about 100 hrs on 16 3.1-GHz Intel Xeon processors**
**Consider simulation** **time interval 1200-7200 seconds**
**RBGen software (POD with Anasazi eigensolver)**
- Compute 151 basis vectors for 151 snapshots (687K length snapshot): 600 sec on 32 3.1-GHz Intel Xeons
**Compare snapshots intervals of 20 and 40 seconds**
- Much cheaper if snapshots over larger intervals provide sufficient information.

## Future Anasazi Development

- New solvers: Generalized Davidson, ???
- Implement StatusTest mechanism:
- Convergence
- Krylov subspace restarting
- Switching to a different solver
- and much, much more!
- Abstract interface to orthogonalization.
- Support for complex scalar types.
- Abstract interface for Rayleigh-Ritz process.
- Solver managers:
- Contains logic not essential to an eigensolver iteration.
- Orchestrates multiple solvers (e.g., hybrid solver).
- Guarantees solution properties.
- Python interface.

## Concluding Remarks(audience participation)

- Are there any additional capabilities that you may need from a eigensolver?
- Do you have any questions, comments, or suggestions about the user interface?
- Are there any other solvers that you would like to see in Anasazi?