# Presentation: 001001

## An Overview of Belosand AnasaziOctober 16th, 11-12amHeidi ThornquistTeri BarthRich LehoucqMike HerouxComputational Mathematics and Algorithms

- 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
- AztecOO and ARPACK
- Solver design requirements
- Belos and Anasazi
- Current Status
- Current Development
- Teuchos
- Arbitrary precision BLAS
- Concluding remarks

## Background

- Several iterative linear solver / eigensolver libraries exist:
- PETSc, SLAP, LINSOL, Aztec(OO), ...
- LOBPCG (hypre), ARPACK, ...
- None are written in C++
- None of the linear solver libraries can efficiently deal with multiple right-hand sides
- Current eigensolver libraries offer one algorithm for solving symmetric/non-symmetric eigenproblems

## AztecOO Linear Solver Library

- 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.
- Uses Epetra to perform underlying vector space operations.
- Interface to matrix-vector product is defined by the user through the EpetraOperator.

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

## Some design requirements for Linear Solvers /Eigensolvers

- Abstract interface to an operator-vector product (Op(A)v) and preconditioning (if necessary).
- Ability to allow the user to enlist any linear algebra package for the elementary vector space operations essential to the algorithm
- Ability to allow the user to define convergence of any algorithm (a.k.a. status testing).
- Ability to allow the user to determine the verbosity level, formatting, and processor for the output.

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

**	****while ( ****myStatusTest.GetStatus(****...****)**
**	****	****	****!= converged )**
**     ...**
**     ...**
**     ****% Compute operator-vector product **
**     ****myProblem.****ApplyOp(****NOTRANS,****v,w****)****; **
**     ******** ****= ****w****.dot(****v****)****;**
**     ...**
**	****end**
**	****myIO.print(****something****)****;**
**  ****	****return (****solution****);**
**end**

**Generic Algorithm**

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

## New Packages:  Belos and Anasazi

- Next generation linear solvers (Belos) and eigensolvers (Anasazi) libraries, written in templated C++.
- Provide a generic interface to a collection of algorithms for solving large-scale linear problems and eigenproblems.
- Algorithm implementation is accomplished through the use of abstract base classes.  Interfaces are derived from these base classes to operator-vector products, status tests, and any arbitrary linear algebra library.
- Includes _block_ linear solvers (GMRES, CG, etc.) and eigensolvers (IRAM, etc.).
- Release:  Spring/Summer 2004

## Why are Block Solvers Useful?(from an applications point of view)

- Block Eigensolvers ( Op(A)X = X ):
- Block Linear Solvers ( Op(A)X = B ):

-   Reliably determine multiple and/or clustered eigenvalues.
-   Example applications:
-   Stability analysis
-   Bifurcation analysis (LOCA)

-   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

## Why are Block Solvers Useful?(from a computational point of view)

- _Remember:_  These are iterative solvers!
- Solutions for block solvers are obtained through operator-multivector computations (Op(A)V).
- However, single-vector iterative solvers use operator-vector computations (Op(A)v).
- _What’s the difference?_

- BLAS Performance
-  Op(A)v is Level 1
-  Op(A)V is Level 2

## Current Status(Trilinos Release 3.2d)

- Belos:
- Solvers:  BlockCG and BlockGMRES
- Interfaces are defined for Epetra and NOX/LOCA
- BelosEpetraOperator allows for integration into other codes
- Epetra interface accepts EpetraOperators, so can be used with ML, AztecOO, Ifpack, Belos, etc...
- Block size is independent of number of right-hand sides
- Anasazi:
- Solver:  Implicitly Restarted Arnoldi
- Interfaces are defined for Epetra and NOX/LOCA
- Can solve standard and generalized eigenproblems
- Epetra interface accepts EpetraOperators, so can be used with Belos, AztecOO, etc...
- Block size is independent of number of requested eigenvalues

## Current Belos/Anasazi Development

- LinearSolver / Eigensolver Class
- LinearProblem / Eigenproblem Class
- StatusTest Class
- IOManager Class
- Incorporate parameter lists ( developed in Teuchos ).
- Replace current Belos/Anasazi abstract linear algebra classes with TSFCore abstract linear algebra classes.

## Current Belos/Anasazi DevelopmentLinearSolver / Eigensolver Class

- Provide an abstract interface to Belos/Anasazi solvers.
- Methods:
- Set Methods:  SetProblem, SetParameters, SetStatusTester
- Get Methods:  GetProblem, GetParameters, GetStatusTester, currIteration
- Solve Methods:  Solve, Iterate
- Auxilliary Method(s):  Reset

## Current Belos/Anasazi DevelopmentLinearProblem / Eigenproblem Class

_**Belos**_
- Allows preconditioning, scaling, etc. to be removed from the algorithm.
- Methods:
- Set Methods:  SetOperator, SetLHS, SetRHS, SetL/RScale, SetL/RPrecond
- Get Methods:  GetOperator, GetLHS, GetRHS, GetL/RPrecond
- Query Methods:  IsL/RPrecond
- Apply Methods: ApplyOperator

_**Anasazi**_
- Differentiates between standard and general eigenvalue problems.
- Allows spectral transformations to be removed from the algorithm.
- Methods:
- Set Methods:  SetInitVec, SetOperator, SetMatrixA, SetMatrixB
- Get Methods:  GetInitVec, GetOperator, GetMatrixA, GetMatrixB, GetEvals, GetEvecs, GetResids
- Apply Methods:  ApplyOperator, ApplyMatrixA, ApplyMatrixB
- IP/Norm Methods:  AInProd, BInProd, BMvNorm

## Current Belos/Anasazi DevelopmentStatusTest Class

- Similar to NOX / AztecOO.
- Allows user to decide when solver is finished.
- Multiple criterion can be logically connected.
- Abstract base class methods:
**StatusType CheckStatus( int currIteration, ****	****	****	****	****	****	****LinearProblem<Scalar>& myProblem );**
**StatusType CheckStatus( int currIteration,****	****	****	****	****	****	****Eigenproblem<Scalar>& myProblem );**
**StatusType GetStatus() const;**
**ostream& Print( ostream& os ) const;**
- Derived classes
- Maximum iterations/restarts, residual norm, combinations {AND,OR}, ...

## Current Belos/Anasazi DevelopmentIOManager Class

- Allows user to control I/O processor, precision, and verbosity.
- Methods:
- Set Methods:
**void reset(Teuchos::ParameterList& myIOParams );**
- Query Methods:
**bool isPrintProc() const;**
**bool isPrintLevel() const;**
**bool isPrintProcLevel() const;**
- Output Methods:  overload “<<“ operator using given precision.

## Future Belos/Anasazi Algorithms

_**Belos**_
- CG/FCG
- GMRES/FGMRES
- BiCGSTAB

_**Anasazi**_
- LOBPCG
- Restarted Preconditioned Eigensolver (RPE)
- Generalized Davidson

## Future Belos/Anasazi Features

_**Belos**_
- Preprocessing capabilities (CG)
- User determined Krylov-subspace breakdown tolerances

_**Anasazi**_
- Preprocessing capabilties
- User determined sorting
- User determined Krylov-subspace breakdown tolerances

## Belos / Anasazi Development(audience participation)

- Are there any capabilities you may need from a linear/eigensolver that have not been addressed?
- Do you have any questions, comments, or suggestions about the proposed user interface?
- Are there any other solvers that you think would be useful?

## Teuchos

**Teuchos**

- Optional, lightweight utility package
- Takes advantage of advanced features of C++:
- Templates
- Standard Template Library (STL)
- Planned features:
- Parameter list of arbitrary types (NOX)
- Memory management package and error-handler  (TSFCoreUtils)
- Templated interfaces to BLAS, LAPACK, dense matrix (TPetra)
- Command line and XML parsers, random number generators, and timers (TSF)

## BLAS

- Basic Linear Algebra Subroutines
- Robust, high-performance dense vector-vector, matrix-vector, and matrix-matrix routines
- Written in Fortran
- Fast and efficient
- Fortran is highly portable; however, the interface between it and C++ is platform-dependent
- Templated C++ wrappers for the BLAS Fortran routines are created through autoconf.
- Two main benefits of using a templated interface
- Template Specialization:  template specialization can be used to support the four native datatypes of existing Fortran BLAS code
- General Implementation:  a generic C++ implementation can be created to handle many arbitrary datatypes

## Arbitrary Datatypes

- Any datatype that defines zero, one, addition, subtraction, and multiplication can use nearly all BLAS functions
- Some require square root or division
- Example:  ARPREC
- arbitrary precision library that uses arrays of 64-bit floating-point numbers to represent high-precision floating-point numbers.
- ARPREC values behave just like any other floating-point datatype, except the maximum working precision (in decimal digits) must be specified before any calculations are done
**	****	****	****	****mp::mp_init(200);**
- Illustrate the use of ARPREC with an example using Hilbert matrices.

## Hilbert Matrices

- A Hilbert matrix ***H******N*** is a square ***N***-by-***N*** matrix such that:
- Notoriously ill-conditioned
**κ****(*****H*****3****)**  524
**κ****(*****H*****5****)**  476610
**κ****(*****H*****10****)**  1.6025 x 1013
**κ****(*****H*****20****)**  7.8413 x 1017
**κ****(*****H*****100****)**  1.7232 x 1020
- Theoretically, Hilbert matrices are SPD

## Cholesky Factorization and Hilbert Matrices

- In practice, floating-point rounding error causes higher-order Hilbert matrices to no longer be positive definite
- Cholesky factorization is often done to determine whether or not a given matrix is positive definite
- As such, a Cholesky factorization will often fail on higher-order Hilbert matrices, even though such matrices “should” be positive definite

| Precision | Largest N for which Cholesky Factorization is successful |
| --- | --- |
| Single Precision | 8 |
| Double Precision | 13 |
| Arbitrary Precision (20) | 29 |
| Arbitrary Precision (40) | 40 |
| Arbitrary Precision (200) | 145 |
| Arbitrary Precision (400) | 233* |

## Cholesky Factorization and Hilbert Matrices

- In addition to providing a higher limit on ***N***, using arbitrary precision reduces error.
- By solving the linear system:                                       where the vector ***b*** is constructed such that the true solution ***x*** is equal to a 1-vector we can look at the forward error.

## Templated BLAS Conclusions

- Adding a templated interface to the BLAS allows the use of arbitrary datatypes
- This is useful in particular for the inclusion of arbitrary-precision arithmetic
- Such arithmetic may be expensive, but the increased precision may be required for ill-conditioned problems
- Arbitrary-precision arithmetic can be used in only the places where it is needed most in order to retain speed
- Templated tools libraries are essential for the development of templated solver libraries.