# Presentation: 153123

## Teuchos:  Utilities for Developers    &    Users  November 2nd, 3:30-4:30pmRoscoe BartlettMike HerouxKris KampshoffKevin LongPaul SextonHeidi Thornquist

**Teuchos:  Utilities ****for**** ****Developers    &    Users****  November 2****nd****, 3:30-4:30pm****Roscoe Bartlett****Mike Heroux****Kris Kampshoff****Kevin Long****Paul Sexton****Heidi Thornquist**

- 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?

## What is Teuchos (TEF-hos)?

- Lightweight, portable utility package
- What’s new?
- First release in Trilinos 4.0
- Templated BLAS/LAPACK interface
- Templated serial dense matrix / vector
- Timing classes
- Templated parameter list
- Numeric traits mechanisms ( ordinal / scalar )
- Memory management:
- Reference-counted pointers
- Pre-allocated workspace
- Exception handling
- Command-line parsing
- “Extended” utilities:
- XML parsing
- MPI communicators / traits

- Users

- Developers

## Templated Parameter List

- Design based on NOX parameter list
- Adopted by Amesos, ML, Ifpack, AztecOO, ...

-    Parameters can be numeric datatypes, pointers to vectors or             functions, and other parameter lists.

## Numeric Traits Mechanisms

- _Assumed_:  addition, subtraction, multiplication, and division
- OrdinalTraits
- zero & one
- int & long int
- ScalarTraits
- zero, one, magnitude type, absolute value, conjugate,
-      square root, random number generator, ...
- std::numeric_limits
- float, double, complex<float>, and complex<double>
- _Required_: none, a compile-time check is performed.
- What about arbitrary precision arithmetic?

- Generic programming technique that uses templated interfaces to define the standard behavior of datatypes.

## Arbitrary Precision LibrariesARPREC & GNU MP

- ARPREC (Bailey, et. al.)
- uses arrays of 64-bit floating-point numbers to represent high-precision floating-point numbers.
- provides support for datatypes with differing precision levels.
**--enable-teuchos-arprec**
- GNU MP
- uses fullwords, arbitrary precision mantissa, and limited precision exponent to represent high-precision floating-point numbers.
- platform specific kernels written in assembly language perform common inner-loop functions.
- provides support for datatypes with differing precision levels.
**--enable-teuchos-gmp**
- Caveat:  Arbitrary precision arithmetic is expensive, but can be used selectively to overcome numerical inaccuracy.

## Templated BLAS Interface

- Basic Linear Algebra Subroutines (BLAS)
- 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.

## Templated BLAS Interface

**Using ARPREC**

**Using Native Datatypes**

**Teuchos::BLAS<int,double> blas;**

**mp::mp_init(200);**
**Teuchos::BLAS<int,mp_real> blas;**

**blas.GEMM(...)**

**blas.GEMV(...)**

**blas.AXPY(...)**

**blas.AXPY(...)**

**blas.GEMV(...)**

**blas.GEMM(...)**

**...**

**...**

**...**

**...**

## Templated LAPACK Interface

- Linear Algebra PACKage (LAPACK)
- Robust, high-performance dense routines for solving systems of simultaneous linear equations, least-squares solutions of linear systems of equations, eigenvalue problems, and singular value problems.
- Templated C++ wrappers for the LAPACK Fortran routines are created for the four native datatypes.
- no general implementations at this time
- Templated serial dense matrix / vector class designed to be compatible with BLAS/LAPACK interfaces.
- Fortran-style storage
- Similar to Epetra_SerialDenseMatrix

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

## Timing Classes

- Teuchos::Time
- Basic wall-clock timing class
- Teuchos::TimeMonitor
- Uses Teuchos::Time internally
- Constructor starts timer, destructor stops it
- Can also maintain an array of timers

**Teuchos::Time timer(“My Timer”);**
**timer.start();**
**...**
**timer.stop();**
**double comp_time = timer.totalElapsedTime()**

## Teuchos Users Overview

- Parameter list is easy to use and quickly being adopted by most Trilinos packages.
- Numeric traits mechanism allows the use of arbitrary datatypes in templated code.
- 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.
- See Teuchos website for more info:
**http://software.sandia.gov/Trilinos/packages/teuchos**