# Presentation: 458971

## An Overview of the Trilinos Package ArchitectureMichael A. HerouxSandia National LaboratoriesTrilinos Workshop at Copper MountainMarch 30, 2004

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

## The Team

**Ross Bartlett **
- Lead Developer of TSFCore
**Paul Boggs **
- Developer of TSF
**Jason Cross **
- Developer of Jpetra
**David Day **
- Developer of Komplex
**Bob Heaphy **
- Lead developer of Trilinos SQA
**Mike Heroux **
- Trilinos Project Leader
- Lead Developer of Epetra, AztecOO,
- Kokkos, Komplex and IFPACK, TSF
- Developer of Amesos, Belos** **
**Robert Hoekstra **
- Developer of Epetra
**Russell Hooper**
- Developer of NOX

**Vicki Howle **
- Developer of Belos and TSF
**Jonathan Hu **
- Developer of ML
**Tammy Kolda **
- Lead Developer of NOX** **
**Rich Lehoucq **
- Developer of Anasazi and Belos
**Kevin Long **
- Lead Developer of TSF,
- Developer of Belos
**Roger Pawlowski **
- Lead Developer of NOX
**Michael Phenow **
- Trilinos Webmaster** **
**Eric Phipps **
- Developer of LOCA and NOX
**Andrew Rothfuss **
- Developer of TSF** **

**Andrew Salinger **
- Lead Developer of LOCA** **
**Marzio Sala**
- Lead author of Trilinos Tutorial
- Developer of ML and Amesos
**Paul Sexton **
- Developer of Epetra and Tpetra
**Ken Stanley **
- Lead Developer of Amesos
**Heidi Thornquist **
- Lead Developer of Anasazi and Belos
**Ray Tuminaro **
- Lead Developer of ML
**Jim Willenbring **
- Developer of Epetra and Kokkos.
- Trilinos library manager** **
**Alan Williams **
- Developer of Epetra

## Outline of Talk

- Motivation
- Trilinos Package Concepts.
- Overview of Major Packages.
- Examples Using Trilinos.
- SQA/SQE.
- Availability and support.
- Concluding remarks.

## Motivation For Trilinos

- Sandia does A LOT of solver work.
- Sandia uses MANY external solver packages.
- When I started at Sandia in May 1998:
- Aztec was a mature package.  Used in many codes.
- FETI, PETSc, DSCPack, Spooles, ARPACK, DASPK, and many other codes were (and are) in use.
- New projects were underway or planned in multi-level preconditioners, eigensolvers, non-linear solvers, etc...
- New application capabilities were underway, demanding new solvers.
- The challenges:
- Little or no coordination was in place to:
- Efficiently reuse existing solver technology.
- Leverage new development across various projects.
- Support solver software processes.
- Provide consistent solver APIs for applications.
- ASCI was forming software quality assurance/engineering (SQA/SQE) requirements:
- Daunting requirements for any single solver effort to address alone.

## Evolving Trilinos Solution

**Evolving Trilinos Solution**

- Trilinos1 is an evolving framework to address these challenges:
- Includes common core set of vector, graph and matrix classes (Epetra).
- Provides a common abstract solver API (TSF).
- Provides a ready-made package infrastructure:
- Source code management (cvs, bonsai).
- Build tools (autotools).
- Automated regression testing (queue directories within repository).
- Communication tools (mailman mail lists).
- Specifies requirements and suggested practices for SQA.
- In general allows us to categorize efforts:
- Efforts best done at the Trilinos level (useful to most or all packages).
- Efforts best done at a package level (peculiar or important to a package).
**Allows package developers to focus only on things that are unique to their package.**

- 1. Trilinos loose translation: “A string of pearls”

## Trilinos Benefits: The Three “I”s

**Infrastructure:**
- Repository (CVS)
- Issue Tracking (Bugzilla)
- Communication (Mailman)
- Debugging (Bonsai)
- Jumpstart (new_package)
- Portable build process (Autotools: configure, make)
- SQA Tools and Policies.
- Documentation: Developers Guides, Installation Guide, Tutorial, etc.
**Interfaces:**
- Interoperability: Between Trilinos Packages, with external SW.
- Extensible.
- Without Interdependence.
**Implementations:**
- Real, working code that implements all interfaces.
- Solid default implementations, but *not required to use*.

## Trilinos Package Concepts

## Trilinos Packages

- Trilinos is a collection of *Packages*.
- Each package is:
- Focused on important, state-of-the-art algorithms in its problem regime.
- Developed by a small team of domain experts.
- Self-contained: No explicit dependencies on any other software packages (with some special exceptions).
- Configurable/buildable/documented on its own.
- Sample packages: NOX, AztecOO, IFPACK, Meros.
- Special package collections: Petra, TSF, Teuchos.

##  Notes:  

-  Notes:
- ASCI Algorithms funds much of Trilinos development (LDRD, CSRF, MICS also).
- All packages available (except TOX, aka Rhythmos).
- All information available at Trilinos website:
- 	_software.sandia.gov/trilinos_

## Slide 10

## Three Special Trilinos Package Collections

**Petra**: Package of concrete linear algebra classes: Operators, matrices, vectors, graphs, etc.
- Provides working, parallel code for basic linear algebra computations.
**TSF**: Packages of abstract solver classes: Solvers, preconditioners, matrices, vectors, etc.
- Provides an application programmer interface (API) to any other package that implements TSF interfaces.
- Inspired by HCL.
**Teuchos (pronounced Tef-hos)**: Package of basic tools:
- Common Parameter list, smart pointer, error handler, timer.
- Interface to BLAS, LAPACK, MPI, XML, ...
- Common traits mechanism.
- Goal: Portable tools that enhance interoperability between packages.

## Dependence vs. Interoperability

- Although most Trilinos packages have no explicit dependence, each package must interact with *some* other packages:
- NOX needs operator, vector and solver objects.
- AztecOO needs preconditioner, matrix, operator and vector objects.
- Interoperability is enabled at configure time.  For example, NOX:
- --enable-nox-lapack     compile NOX lapack interface libraries
- --enable-nox-epetra     compile NOX epetra interface libraries
- --enable-nox-petsc      compile NOX petsc interface libraries
- Trilinos is a vehicle for:
- Establishing interoperability of Trilinos components...
- Without compromising individual package autonomy.
- Trilinos offers five basic interoperability mechanisms.

## Trilinos Interoperability Mechanisms

- M1: *Package* accepts user data as Epetra or TSF objects.
- =>Applications using Epetra/TSF can use *package*.
- M2: *Package* can be used via TSF abstract solver classes.
- => Applications or other packages using TSF can use *package*.
- M3: *Package* can use Epetra for private data.
- =>* Package * can then use other packages that understand Epetra.
- M4: *Package* accesses solver services via TSF interfaces.
- =>* Package * can then use other packages that implement TSF interfaces.
- M5: *Package* builds under Trilinos configure scripts.
- =>* Package* can be built as part of a suite of packages.
- => Cross-package dependencies can be handled automatically.

##   Interoperability Example: AztecOO

- AztecOO: Preconditioned Krylov Solver Package.
- Primary Developer: Mike Heroux.
- Minimal *explicit, essential* dependence on other Trilinos packages.
- Uses abstract interfaces to matrix/operator objects.
- Has independent configure/build process (but can be invoked at Trilinos level).
- Sole dependence is on Epetra (but easy to work around).
*Interoperable* with other Trilinos packages:
- Accepts user data as Epetra matrices/vectors.
- Can use Epetra for internal matrices/vectors.
- Can be used via TSF abstract interfaces.
- Can be built via Trilinos configure/build process.
- Can provide solver services for NOX.
- Can use IFPACK, ML or AztecOO objects as preconditioners.

## Trilinos Package Categories

- TSF

- Epetra

- NOX

- Nonlinear solvers

- Vector, graph, matrix service classes

- Abstract solver API

- Preconditioned Krylov solvers

- Algebraic Preconditioners

- Multi-level Preconditioners

## Trilinos Package Interoperability

## What Trilinos is not

- Trilinos is not a single monolithic piece of software. Each package:
- Can be built independent of Trilinos.
- Has its own self-contained CVS structure.
- Has its own Bugzilla product and mail lists.
- Development team is free to make its own decisions about algorithms, coding style, release contents, testing process, etc.
- Trilinos top layer is not a large amount of source code:
- Trilinos repository contains 452,187 source lines of code (SLOC).
- Sum of the packages SLOC counts :	445,937.
- Trilinos top layer SLOC count:	    	    6, 250 (1.4%).
- Trilinos is not “indivisible”:
- You don’t need all of Trilinos to get things done.
- Any collection of packages can be combined and distributed.
- Current public release contains only 6 of the 20+ Trilinos packages.

## Overview of Trilinos Packages

## Trilinos Concrete Support Component: Petra

- 1Petra is Greek for “foundation”.

- Petra1 provides distributed matrix and vector services.
- Exists in basic form as an object model:
- Describes basic user and support classes in UML, independent of language/implementation.
- Describes objects and relationships to build and use matrices, vectors and graphs.
- Has 3 implementations under development.

## Petra Implementations

- Three version under development:
- Epetra (*Essential Petra*):
- Current production version.
- Restricted to real, double precision arithmetic.
- Uses stable core subset of C++ (circa 2000).
- Interfaces accessible to C and Fortran users.
- Tpetra (*Templated Petra*):
- Next generation C++ version.
- Templated scalar and ordinal fields.
- Uses namespaces, and STL: Improved usability/efficiency.
- Jpetra (*Java Petra*):
- Pure Java. Portable to any JVM.
- Interfaces to Java versions of MPI, LAPACK and BLAS via interfaces.

## Epetra

- Package of _concrete linear algebra classes_:
- 		Operators, matrices, vectors, graphs, etc.
- Working, parallel code for basic linear algebra computations.
- Uses stable core subset of C++
- C/Fortran wrappers
- Restricted to real, double precision arithmetic
- Concrete implementation of the Petra object model

## Epetra User Class Categories

- Sparse Matrices: 	***RowMatrix***, (CrsMatrix, VbrMatrix, FECrsMatrix, FEVbrMatrix)
- Linear Operator:	***Operator*****: **(AztecOO, ML, Ifpack)
- Dense Matrices:	DenseMatrix, DenseVector, BLAS, LAPACK,			SerialDenseSolver
- Vectors: 		Vector, MultiVector
- Graphs:		CrsGraph
- Data Layout: 		Map, BlockMap, LocalMap
- Redistribution:	Import, Export, LbGraph, LbMatrix
- Aggregates: 	 	LinearProblem
- Parallel Machine: 	***Comm***, (SerialComm, MpiComm, MpiSmpComm)
- Utilities: 	   	Time, Flops

## Summary of Epetra Features

- Basic Stuff: What you would expect.
- Variable block matrix data structures.
- Multivectors.
- Arbitrary index labeling.
- Flexible, versatile parallel data redistribution.
- Language support for inheritance, polymorphism and extensions.
- View vs. Copy.

## Trilinos Solver Framework (TSF)

- Epetra, AztecOO, Ifpack, ML, etc.PETSc, SuperLU, Hypre, HSL,ScaLapack
- TSF is an abstract class hierarchy:
- Inspired by HCL, especially VectorSpace concept.
- Provides uniform API to solvers, vectors, matrices.
- Allows integration of many solvers via implementation of abstract classes.
- Supports “generic” programming.
- Provides compositional classes.
- Composed of TSFExtended, TSFCore, TSFCoreUtils.

- Lots of good solvers available

## Generic Programming using TSF

- Generic Programming: 	Implementation of algorithms using abstract 	interfaces.
- Example: CG using TSF interfaces.
- Allows use of CG with *any* vector/matrix classes that implement TSF interfaces.
- Very powerful for complex algorithms: Robust Block GMRES, etc.

## Aggregate Operator Construction

- TSF facilitates implicit (and explicit) construction of operators:
- Partitioned (block):
- Composite:
- Sum:
- Inverse:
- Others: Zero, Identity, Transpose, ...
- Recursively.

## AztecOO

- Aztec is the workhorse solver at Sandia:
- Extracted from the MPSalsa reacting flow code.
- Installed in dozens of Sandia apps.
- 1700+ external licenses.
- AztecOO leverages the investment in Aztec:
- Uses Aztec iterative methods and preconditioners.
- AztecOO improves on Aztec by:
- Using Epetra objects for defining matrix and RHS.
- Providing more preconditioners/scalings.
- Using C++ class design to enable more sophisticated use.
- AztecOO interfaces allows:
- Continued use of Aztec for functionality.
- Introduction of new solver capabilities outside of Aztec.

## ML: Multi-level Preconditioners

- ML package developed by Ray Tuminaro, Jonathan Hu, Marzio Sala.
- Critical technology for scalable performance of some key apps.
- ML compatible with other Trilinos packages:
- Accepts user data as Epetra_RowMatrix object (abstract interface).
- Any implementation of Epetra_RowMatrix works.
- Implements the Epetra_Operator interface.
- Allows ML preconditioners to be used with AztecOO and TSF.
- Can also be used completely independent of other Trilinos packages.

## IFPACK: Algebraic Preconditioners

- Overlapping Schwarz preconditioners.
- Accept user matrix via abstract matrix interface (Epetra versions).
- Uses Epetra for basic matrix/vector calculations.
- Supports simple perturbation stabilizations and condition estimation.
- Separates graph construction from factorization, improves performance substantially.
- Compatible with AztecOO and TSF. Can be used with NOX.

## Komplex: Complex linear solver

- Most algorithms work for complex numbers (with real numbers as a special case).
- Majority of our applications produce real-valued data.
- Solver development has been focused on real-valued problems.
- Writing complex versions of all software is not appealing.
- Alternative:  Consider equivalent real formulations (ERFs).
- Komplex is an add-on module to AztecOO that:
- Intelligently builds an ERF for a complex valued problem.
- Computes the real-valued solution using AztecOO.
- Returns the complex result to user.
- This is an effective approach in important practical settings.

- (A+*i*B)(x+*i*y) = (b+*i*c)

- A  -B    x      b
-       =
- B  A     y      c

## Komplex Formulation

- Consider a complex-valued matrix *C*:
- With each entry:

**Rewrite as real-valued of twice the dimension:**

## NOX: Nonlinear Solvers

- Suite of nonlinear solution methods:
- Uses abstract vector and “group” interfaces.
- Allows flexible selection and tuning of various strategies:
- Directions.
- Line searches.
- Epetra/AztecOO/ML, LAPACK, PETSc implementations of abstract vector/group interfaces.
- Designed to be easily integrated into existing applications.

## LOCA

- Library of continuation algorithms.
- Provides
- Zero order continuation
- First order continuation
- Arc length continuation
- Multi-parameter continuation (via Henderson's MF Library)
- Turning point continuation
- Pitchfork bifurcation continuation
- Hopf bifurcation continuation
- Phase transition continuation
- Eigenvalue approximation (via ARPACK or Anasazi)

## Amesos: Direct Solver Wrappers

- Direct sparse solver use at Sandia:
- Salinas (Structures): DSCPACK (Raghavan), SPOOLES (Ashcraft), others.
- Xyce (Circuits): Kundert, SuperLU serial.
- PCx (LP): DSCPACK, WSMP.
- Numerous other uses.
- Amesos contains wrapper classes to important third party direct sparse solvers:
- Use Epetra objects.
- Provide data redistribution capabilities (e.g., replication).
- Provide common look-and-feel across variety of solvers.
- Provide common resource for direct solvers at Sandia.
- Presently wraps: DscPack, KLU, Mumps, SuperLU, UMFPack, Kundert
- Goal:  “Make solving *Ax=b* as easy as it sounds.”

## Epetraext: Extensions to Epetra

- Library of useful classes not needed by everyone.
- Most classes are types of “transforms”.
- Examples:
- Graph/matrix view extraction.
- Epetra/Zoltan interface.
- Explicit sparse transpose.
- Singleton removal filter.
- Static condensation filter.
- Overlapped graph constructor.
- Graph colorings.
- Permutations.
- Sparse matrix-matrix multiply.
- Matlab, MatrixMarket I/O functions.
- ...
- Most classes are small, useful, but non-trivial to write.

## Meros

- Meros: Preconditioner package for incompressible NS problems.
- Addresses problems: *Ax = b.*
- where:
- Makes use of TSF to orchestrate use of:
- ML, Epetra, AztecOO, Ifpack.
- Provides rapidly-developed, scalable implementation of state-of-the-art preconditioner.
- Results shown below.

## 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 matrix-vector products, status tests, and any arbitrary linear algebra library.

## New Package: Kokkos

- Goal:
- Isolate key non-BLAS kernels for the purposes of optimization.
- Kernels:
- Dense vector/multivector updates and collective ops (not in BLAS).
- Sparse MV, MM, SV, SM.
- Serial-only for now.
- Reference implementation provided.
- Mechanism for improving performance:
- Default is aggressive compilation of reference source (Fortran loops).
- BeBOP: Jim Demmel, Kathy Yelick, Rich Vuduc, UC Berkeley.
- Vector version: Cray.

## Teuchos

- Utility package of commonly useful tools:
- ParameterList class: key/value pair database, recursive capabilities.
- LAPACK, BLAS wrappers (templated on ordinal and scalar type).
- Dense matrix and vector classes (compatible with BLAS/LAPACK).
- FLOP counters, Timers.
- Ordinal, Scalar Traits support: Definition of ‘zero’, ‘one’, etc.
- Reference counted pointers, and more...
- Takes advantage of advanced features of C++:
- Templates
- Standard Template Library (STL)
- ParameterList:
- Allows easy control of solver parameters.

**Notes:**
Teuchos is being designed as an optional utility package that takes advantage of some of the more advanced features of C++, including templates and the STL.
Planned features for Teuchos include: a parameter list of arbitrary types, a memory management package, and templated interfaces to the BLAS, LAPACK and a dense matrix class.
My talk will focus on the templated interface to the BLAS. First, I’ll give you a brief background on the BLAS and C++ templates.

## Teuchos BLAS/LAPACK

- Teuchos provides C++ wrappers for BLAS/LAPACK Fortran routines
- Insulates users of Teuchos from details of Fortran function calls
- Instead of using simple wrappers, Teuchos provides templated wrappers
- Two main benefits of using a templated interface
- Template Specialization
- Fortran BLAS supports 4 datatypes:
- Single- and double- precision real, single- and double-precision complex
- We can use template specialization to support these datatypes using existing Fortran BLAS code
- General Implementation
- We can write a generic implementation in C++ that can handle many arbitrary datatypes

**Notes:**
Teuchos provides what are called wrappers for the BLAS Fortran routines. The wrapper handles all translation between C++ and Fortran and insulates the rest of Teuchos from the details.
Using templated wrappers adds another layer of abstraction, and this fits naturally with the BLAS. You’ll see why in a moment after I have given a brief background on templates.

## NewPackage Package

- NewPackage provides jump start to develop/integrate a new package.
- NewPackage is a “Hello World” program and website:
- Simple but it does work with autotools.
- Compiles and builds.
- NewPackage directory contains:
- Commonly used directory structure: src, test, doc, example, config.
- Working autotools files.
- Documentation templates (doxygen).
- Working regression test setup.
- Substantially cuts down on:
- Time to integrate new package.
- Variation in package integration details.
- Development of website.

## What do algorithms developers get from the Trilinos package architecture?

## Pliris

- Migration of Linpack benchmark code to supported environment.
- Used NewPackage as starting point.
- Put under source control for the first time.
- Portable build process.
- Can accept linear problem as Epetra data.
- Has its own Bugzilla product, mail lists.
- Source code browsable via Bonsai.
- Portable interface to BLAS/LAPACK.
- And more...
- But Pliris is still an independent piece of code...and better off after the process.

## Trilinos Package Interoperability (Release 4)

- Based on this chart:
-  AztecOO depends on Epetra, but Epetra is independent of AztecOO
-  NOX can use Epetra, but is independent of Epetra.

| 
Package |  | Depends On |  |  |  | Can  Use |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  | Epetra | Epetra Ext | AztecOO | Komplex | Ifpack | Amesos | ML | NOX | TSF Core | TSF Ext | Belos | Meros | Anasazi |
| Epetra |  |  |  |  |  |  |  |  |  |  |  |  |  |
| Epetra Ext |  |  |  |  |  |  |  |  |  |  |  |  |  |
| AztecOO |  |  |  |  |  |  |  |  |  |  |  |  |  |
| Komplex |  |  |  |  |  |  |  |  |  |  |  |  |  |
| Ifpack |  |  |  |  |  |  |  |  |  |  |  |  |  |
| Amesos |  |  |  |  |  |  |  |  |  |  |  |  |  |
| ML |  |  |  |  |  |  |  |  |  |  |  |  |  |
| NOX |  |  |  |  |  |  |  |  |  |  |  |  |  |
| TSFCore |  |  |  |  |  |  |  |  |  |  |  |  |  |
| TSF Ext |  |  |  |  |  |  |  |  |  |  |  |  |  |
| Belos |  |  |  |  |  |  |  |  |  |  |  |  |  |
| Meros |  |  |  |  |  |  |  |  |  |  |  |  |  |
| Anasazi |  |  |  |  |  |  |  |  |  |  |  |  |  |

## Applications

## Block preconditioners for incompressible flow

where

Easier to apply MG on blocks than whole system
					  scalability!

Close to traditional pressure-Poisson type solvers 
Approximation choices  different methods

**Block preconditioners for incompressible flow**
- where
- _Easier to apply MG_ on blocks than whole system
- 					  _**scalability!**_
- Close to traditional pressure-Poisson type solvers
- Approximation choices  different methods

**Meros**

## 3D Driven Cavity Results

**3D Driven Cavity Results**

| Grids | 163 | 323 | 643 |
| --- | --- | --- | --- |
|  its @ Re=100 | 15 | 17 | 18 |

| Re | CFL | Oseen
(K&L) | Newton (Oseen K&L) | Newton
(BJB) |
| --- | --- | --- | --- | --- |
| 500 | 5000 | 5 | 24 | 4 |
|  | 10000 | 5 | 36 | 6 |

_**IMPACT**_** Block preconditioning is a practical & efficient way to achieve scalable linear solvers.**
_**Potential Applications**_**:  ****Aria, Charon, Fuego, Goma, MPSalsa, Sundance**

- CFL insensitivity
- h-independence

**3D Driven Cavity**

- Re sensitivity?

**Notes:**
N=64, N3 cells: approximately 4N3 dof.
 NS = ave. # nonlinear steps per time step
 outer = ave. # outer linear iterations per NS
 F = ave. # F iters per outer solve (ave of 16 Ap iters per outer solve)
N=64, N3 cells: approximately 4N3 dof.
 NS = ave. # nonlinear steps per time step
 outer = ave. # outer linear iterations per NS
 F = ave. # F iters per outer solve (ave of 16 Ap iters per outer solve)

From Vicki Howle:

I'd say that the main point is the h-independence. That is the big gain
and the block preconditioners help us get that because applying MG on the
individual blocks is easier then applying it on the whole saddle point
matrix.

The point behind the table you call "CFL insensitivity" is that these
preconditioners are still good for high CFL numbers (i.e., large delta-t).
In particular, we looked at some higher Re problems for which the
steady-state problem was very difficult to solve. However, transient
problems are easier (in the linear algebra), so we could use a transient
formulation of the problem with very large time steps (high CFL) to
converge to the steady solution (a pseudo-transient solver).

The brief answer on "what is BJB" is that it is the same style of block
preconditioner as Kay and Loghin, but with the approximation to the Schur
complement inverse being X^{-1} = -A_p^{-1} (B J B^T) A_p^{-1}. A_p is a
pressure Laplacian term.

Here are some more details on the three columns of that table in case you
want them:

In general, the block preconditioners we are looking at are of the form

  P^{-1} = | F^{-1} 0 | | I  B^T | | I   0     |
           |  0     I | | 0   I  | | 0  X^{-1} |

where X^{-1} is an approximation to the Schur complement inverse S^{-1} =
(-B F^{-1} B^T)^{-1}, and where the original saddle matrix was

 | F  B^T |
 | B   0  |

(We are currently looking at what we need to do in the preconditioning
when there is a stabilization term ----> [ F  B^T; B  C ].)

In the "CFL insensitivity" table, the "Oseen (K&L)" column is showing
results for using the Kay & Loghin preconditioner on the saddle linear
system you get from a discretization of the Oseen equations (i.e., a
Picard iteration on the NS eqns). The Kay and Loghin preconditioner is the
P^{-1} above with the choice of Schur complement approximation X^{-1} =
-F_p A_p^{-1}. The A_p term is a pressure Laplacian (wonderful for MG).
The F_p term is the convection-diffusion operator applied on the pressure
space. Applying the Kay & Loghin preconditioner therefore requires just
one solve on the pressure Laplacian term A_p and one solve on the
convection-diffusion operator F. The draw-back of this method is that the
application code has to provide us with the non-standard operator F_p. (We
are also currently researching ideas for automatic generation of the F_p
term.)

The "Newton (K&L)" column is examining using the Oseen-based Kay & Loghin
preconditioner on the linear system you get from using Newton for your
nonlinear solver. Now the saddle matrix looks like

 | J  B^T |
 | B   0  |

where J is the Jacobian. And the preconditioner is

  P^{-1} = | J^{-1} 0 | | I  B^T | | I   0     |
           |  0     I | | 0   I  | | 0  X^{-1} |

where X^{-1} is still -F_p A_p{-1}. F_p is still the Oseen
convection-diffusion operator on the pressure space.

The idea in this column is to see how this (cheaper) Oseen-based K&L
preconditioner does on the Newton saddle problem. It takes more iterations
since X^{-1} isn't as good an approximation to the inverse of the Newton
Schur complement S = -B J^{-1} B^T.

The final column "Newton (BJB)" is again using Newton for the nonlinear
solver, so our saddle matrix is again

 | J  B^T |
 | B   0  |

In this column, we are using Howard Elman's block preconditioner (which
we've been calling "BJB"). The BJB block preconditioner is of the same
block form as K&L,

  P^{-1} = | J^{-1} 0 | | I  B^T | | I   0     |
           |  0     I | | 0   I  | | 0  X^{-1} |

but with a different choice of Schur complement approximation. For BJB we
have X^{-1} = -A_p^{-1} (B J B^T) A_p^{-1}. Applying this preconditioner
now requires one solve on J and two solves on the pressure Laplacian A_p,
but this X is a better approx to the Newton Schur complement
S = -B J^{-1} B^T. (The arguments behind the derivation of the K&L choice
of X work better on the Oseen problem, so it's not as good a
preconditioner for Newton.) In addition, with the BJB preconditioner
we don't need the non-standard operator F_p.

Hope that makes sense! Let me know if you have any questions.

-vicki

## Inter-operation of Trilinos packages

**Inter-operation of Trilinos packages**

**Nonlinear**
**Solver**

**Linear**
**Solver**

**block**
**precondition**

**NOX****(Epetra)**

**AztecOO**
**(TSF, Epetra)**

**Meros**
**(TSF)**

**MPSalsa**

**(Epetra)**

**Newton-Krylov**
**Methods**

**GMRESR**

**Solution ****component**

**Example ****methods**

**Packages**

**Subblocks**
**(AztecOO, ****ML, Epetra)**

**F****-1****:**
**  ****GMRES/AMG**
**X****-1****:**
**  ****CG/AMG**

**MPSalsa ******** Epetra submatrices**
**Epetra ******** TSF**
**Epetra/TSF matrix add &   ****	****multiply**
**M****-1**** for TSF subblocks****	**

**Notes:**
Say or notate which parts we will implement in TSF (for Salsa problem for example)
Pressures on cell centers, velocities on cell faces
MAC discretization
Implemented in Trilinos/Petra; using ML for multigrid
Outer linear solve: GMRES
F / J solves: GMRES/MG
Ap solves: CG/MG

## Circuit SolversBringing Multiple Tools Together

**Circuit Solvers****Bringing Multiple Tools Together**

- Reduction to “Essential” System
- Constraints (‘Singleton’ Rows)
- Auxiliary Equations (‘Singleton’ Columns)
- Based on Achim Basermann’s Results (NEC Europe)
- Available in EpetraExt.

- Load Balance
- Circuit

- Scale

- RCM

- Singleton
- Filter

- Load Balance
- LinSys

| N | Total Cuts | CondEst | GMRES Iters | LinSolve Time | Newton Steps |
| --- | --- | --- | --- | --- | --- |
| 1220 | ~1000 | 3.00E+05 | 500 | 4.7 | 72 |
| 1220 | 54 | 2.00E+05 | 200 | 1.1 | 45 |
| 1054 | 68 | 1.00E+04 | 127 | 0.36(0.43) | 53 |

- LBC

- LBC+SF+LBL+RCM+SCALE

**Notes:**
Mentioned last year
Now reality in Epetra (M. Heroux)
Example of Composite Transformations in use
Showing reductions in iters, time, steps
Time even with overhead

## ARPREC

- The ARPREC library 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.

**Notes:**
ARPREC is just one such datatype that fits that description.
The ARPREC library uses arrays of floating-point numbers to represent larger floating-point numbers – basically, arrays just get chained together, along with some additional information.
ARPREC values behave just like any other floating-point datatype – that is, it satisfies all the requirements I described in the previous slide – with one caveat: before we use ARPREC values, we need to specify the maximum working precision.

## Hilbert Matrices

- A Hilbert matrix ***H******N*** is a square ***N***-by-***N*** matrix such that:
- For Example:

**Notes:**
Now that I’ve described how the templated BLAS interface works and an arbitrary-precision datatype that we can plug into it, I’ll show how we can take advantage of the increased precision provided by ARPREC. This will require a little more background, so please bear with me; I will eventually get to a point.
A Hilbert matrix is a square matrix where each element is equal to one over the element’s row index plus its column index minus one. For example, this is a three-by-three Hilbert matrix.

## Hilbert Matrices

- 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
- Hilbert matrices introduce large amounts of error

**Notes:**
Hilbert matrices are notoriously ill-conditioned, which means that they introduce large amounts of error into calculations they are involved with. As you can see, as the Hilbert matrices get larger and larger, their condition numbers grow at a phenomenal rate.
The increased precision provided by ARPREC should alleviate much of this error.

## Hilbert Matrices and Cholesky Factorization

- With double-precision arithmetic, Cholesky factorization will fail for ***H******N*** for all ***N*** > 13.
- Can we improve on this using arbitrary-precision floating-point numbers?

| 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+ (out of memory to go higher) |

**Notes:**
For example, with double-precision numbers, Cholesky factorization will fail for any Hilbert matrix with a size greater than 13.
Now, the question is, can we improve on this using arbitrary precision?
The answer is yes. Obviously, the more precision we use, the larger Hilbert matrices we are able to factor. At 200 digits of precision, we are able to factor a 145-by-145 matrix; at 400 digits of precision, we can do even better. This 233 is actually a lower bound, because the machine I was running my tests on actually ran out of memory at this point, so I believe it could probably do a bit better than 233.

## Preliminary Performance Analysis

| N | double | mp_real(40) | Ratio mp_real/double |
| --- | --- | --- | --- |
| 100 | 0.0001257 sec | 0.0437751sec
 | 349 |
| 500 | 0.0045598 sec
 | 1.0673 sec
 | 234 |
| 1000 | 0.093067 sec
 | 4.31996 sec
 | 46.5 |

- Time (sec) of GEMV

## Potential Benefits of Templated Types

- Templated scalar types have great potential:
- Efficient:  Allow expression of algorithms over any abstract field.
- Can facilitate variety of algorithmic studies.
- Allow application developers to study asymptotic behavior of discretizations.
- Can facilitate debugging:  Reduces FP error as source error.
- Use your imagination...

## SQA/SQE

- Software Quality Assurance/Engineering is important.
- Not sufficient to convince ourselves that “we do a good job”.
- At start of development: Must state the practices we will follow.
- At end of development: Must provide evidence that we followed practices.
- Trilinos facilitates SQA/SQE development/processes for packages:
- 32 of 47 ASCI SQE practices are directly handled by Trilinos (no requirements on packages).
- Trilinos provides significant support for the remaining 15.
- Trilinos Dev Guide Part II: Specific to ASCI requirements.
- Trilinos software engineering policies provide a ready-made infrastructure for new packages.
- Trilinos philosophy: Few *requirements*.  Instead mostly *suggested practices.  *Provides package with option to provide alternate process.

## Trilinos Availability/Support

- Trilinos and related packages are available via LGPL.
- Current release (3.1) is “click release”.  Unlimited availability.
- Release 4 scheduled for May 2004.
- New platform facilitates development and support:
- http://software.sandia.gov
- Location of cvs repository, bugzilla, bonzai and mailman servers.
- Accessible from anywhere via ssh/scp.
- Documentation (generated via doxygen) is all available online.

## Mailman Mail Lists

- Each Trilinos package, including Trilinos itself, has four mail lists:
- package-checkins@software.sandia.gov
- CVS commit emails.
- package-developers@software.sandia.gov
- Mailing list for developers.
- package-users@software.sandia.gov
- Issues for package users.
- package-announce@software.sandia.gov
- Releases and other announcements specific to the package.
- Additional list: Trilinos-Leaders@software.sandia.gov
- http://software.sandia.gov/mailman/listinfo/

## Conclusions

- Trilinos services to developers and users:
- The 3 I’s: Infrastructure, Interfaces, Implementations.
- Simplifies installation, support for users of total collection.
- Epetra & TSF promote common APIs across all other Trilinos packages.
- Each package can be built, used independently, and exists as independent project.
**Primary goals: **
**Rapid development and installation of robust numerical solvers.**
**High-quality production software for the critical path.**
- Can SW infrastructure make our jobs easier? Yes!
- NewPackage, mail lists, bugzilla, CVS, Bonsai
- Make you an instant SW engineering professional.
- Package architecture always respects your independence.
- Compatability from augmentation not conformance.
- You are free to “walk away” at any time.

## More information

- http://software.sandia.gov
- http://software.sandia.gov/trilinos
- Additional documentation at my website:http://www.cs.sandia.gov/~mheroux.