# Presentation: 417829

## Scalable Linear Algebra Capability AreaMichael A. HerouxSandia National Laboratories

- 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:**
But just for Tpetra.

## Current Efforts

- Support for multi-precision:
- Double-precision is not always questioned in scientific computing.
- Single-prec. floating point arithmetic can be significantly faster.
- Smaller word size puts less strain on taxed memory hierarchy.
- Multi-precision algorithms allow combination of fast arithmetic and need for higher accuracy.
- Support for newer architectures:
- GPU, CBE, ?
- Can achieve these via a general purpose programming environment with runtime support for desired platforms.
- e.g., Sequoia, RapidMind
- Too much trouble for me.
- Instead, narrow scope to our libraries (i.e., those pesky solvers).

## Tpetra Abstract Interfaces

- We propose a set of abstract interfaces for achieving these goals in the Tpetra library of linear algebra primitives.
- Tpetra is a templated implementation of the Petra Object Model:
- these classes provide data services for many other packages in the Trilinos project (*e.g., linear solvers, eigensolvers, non-linear solvers, preconditioners*)
- successor to Trilinos’s Epetra package
- Tpetra centered around the following interfaces:

-  Comm for providing communication between nodes
-  Map for describing layout of distributed objects.

-  DistObject for redistributing distributed objects
-  Linear algebra object interfaces (Operator, Vector)

## Satisfying Our Goals: Templates

- How do we support multiple data types?
- C++ templating of the scalar type and ordinals.
- Not new, not difficult.
- Compiler support is good enough now.
- This provides generic programming capability, independent of data types.
- Templating implements compile time polymorphism.
- Pro: No runtime penalty.
- Con: Potentially large compile-time penalty.
- This is okay. Compiling is a good use of multicore! :)
- Techniques exist for alleviating this for common and user data types (explicit instantiation)

## Example

**Standard method prototype for apply matrix-vector multiply:**
- template <typename OT, typename ST>
- CrsMatrix::apply(const MultiVector<OT, ST> &x,
-                        MultiVector<OT, ST> &y)

**Mixed precision method prototype (DP vectors, SP matrix):**
- template <typename OT, typename ST>
- CrsMatrix::apply(const MultiVector<OT,ScalarTraits<ST>::dp> &x,
- 		        MultiVector<OT,ScalarTraits<ST>::dp> &y)
**Exploits traits class for scalar types:**
- typename ScalarTraits<ST>::dp;      // double precision w.r.t. ST
- typename ScalarTraits<ST>::hp;      //   half precision w.r.t. ST
- ST ScalarTraits<ST>::one();         // multiplicative identity

**Sample usage in a mixed precision algorithm:**
- Tpetra::MultiVector<int, float> x, y;
- Tpetra::CisMatrix<int, double> A;
- A.apply(x, y);  // SP matrix applied to DP multivector

## C++ Templates

- Example was for float/double but works for:
- complex<float> or complex<double>
- Arbitrary precision (*e.g., GMP, ARPREC*)
- The only requirement is a valid specialization of the traits class.

## The Rest: C++ Virtual Functions

- How do we address our desire to support multiple implementations for these objects?
- C++ virtual functions and inheritance.
- This provides runtime polymorphism.
- Use abstract base classes to encapsulate data and behavior.
- Specific concrete implementations of these interfaces provide adapters to target architectures.
- We will “abstract away” communication, data allocation/placement and computation.

## Tpetra Communication Interface

- Teuchos::Comm is a pure virtual class:
- Has no executable code, interfaces only.
- Encapsulates behavior and attributes of the parallel machine.
- Defines interfaces for basic comm. services between “nodes”, e..g.:
- collective communications
- gather/scatter capabilities
- Allows multiple parallel machine implementations.
- Generalizes Epetra_Comm.
- Implementation details of parallel machine confined to Comm subclasses.
- In particular, Tpetra (and rest of Trilinos) has no dependence on any particular API (*e.g., MPI*).

## Comm Methods

**getRank**()
**getSize**()
- _[barrier](file://Heroux00/u0/mheroux/Trilinos/doc/epetra/html/class_Epetra_Comm.html#a1)_()
**broadcast<Packet>**(Packet *MyVals, int count, int Root)
**gatherAll<Packet>**(Packet *MyVals, Packet *AllVals, int count)
**reduceAll<Packet>**(ReductionOp op,
- 		     int count, const Packet *local, Packet *global)
**scan<Packet>**(ReductionOp op,
- 		int count, const Packet *send, Packet *scans)

**Comm Implementations**

- SerialComm simultaneous supports of serial and parallel coding.
- MpiComm is a thin wrapper around MPI communication routines.
- MpiSmpComm allows use of shared-memory nodes.

## Abstract Node Class

- Trilinos/Kokkos: Trilinos compute node package.
- Abstraction definition in progress.
- Node currently envisioned as an abstract factory class for computational objects.

- Kokkos::Node

- Kokkos::SerialNode

- Kokkos::CudaNode

- Kokkos::TbbNode

- ...

**Example:**
- Kokkos::LocalCrsMatrix<int,double> lclA;
- lclA = myNode.createCrsMatrix(...);
- lclA.submitEntries(...); *// fill the matrix*
- Kokkos::LocalMultiVector<int,double> lclX = myNode.createMV(...),
-                                      lclY = myNode.createMV(...);
- lclA.apply(lclX,lclY); *// apply the matrix operator*

## Abstract Node Class (2)

- Node handles platform-specific details, such as:
- how to allocate memory for the necessary data structures?
- significant in the case of attached accelerators with distinct memory space.
- How to perform the necessary computations?
- Tpetra is responsible only for invoking global communication, via the abstract Comm class.
- In addition to supporting multiple architectures, Tpetra/Kokkos becomes a test bench for research into primitives.
- These abstractions *(hopefully)* provide us with flexibility to tackle a number of platforms.
- Cons:
- m kernels, p platforms  m*p implementations
- Heros can’t improve code they don’t have access to.

## Sample Code Comparison: MV::dot()

**MPI-only****:**
- double dot(int len,
-            double *x,
-            double *y)
- {
-   double lcl = 0.0, gbl;
-   for (int i=0; i<len; ++i)
-     lcl += x[i]*y[i];
-   MPI_ALLREDUCE(lcl,gbl,...);
-   return gbl;
- }

**Tpetra/Kokkos:**
- template <typename ST>
- ST Tpetra::MultiVector<ST>::dot(
-        Comm comm,
-        Kokkos::LocalVector<ST> x,
-        Kokkos::LocalVector<ST> y)
- {
-   Scalar lcl, gbl;
-   lcl = x.dot(y);
-   reduceAll<ST>(comm,SUM,lcl,gbl);
-   return gbl;
- }

- For appropriate choices of Node and Comm, both implementations are equivalent.
- Right hand example is limited only by the available implementations of these classes:
- can determine whether library was compiled with support for GPU, MPI, etc.
- can compose different nodes for heterogeneous

## Example: MPI-Only Ax = b

- App
- Rank 0

- App
- Rank 1

- App
- Rank 2

- App
- Rank 3

- Lib
- Rank 0

- Lib
- Rank 1

- Lib
- Rank 2

- Lib
- Rank 3

- Mem
- Rank 0

- Mem
- Rank 1

- Mem
- Rank 2

- Mem
- Rank 3

- All ranks store *A, x, b* data in locally-visible address space

- Library solves *Ax=b* using distributed memory
- algorithm, communicating via MPI

- App passes matrix and vector values to library data
- classes

- Serial
- Kernel
- Rank 0

- Serial
- Kernel
- Rank 1

- Serial
- Kernel
- Rank 2

- Serial
- Kernel
- Rank 3

- Using: SerialNode, MpiComm

## Example: MPI+Threads Ax = b

**Example: MPI+Threads *****Ax****** ******=****** ******b***

- App
- Rank 0

- App
- Rank 1

- App
- Rank 2

- App
- Rank 3

- Lib
- Rank 0

- Lib
- Rank 1

- Lib
- Rank 2

- Lib
- Rank 3

- Mem
- Rank 0

- Mem
- Rank 1

- Mem
- Rank 2

- Mem
- Rank 3

- Multicore: “PNAS” Layout

- Lib
- Rank 0
- Thread 0   Thread 1   Thread 2  Thread 3

- All ranks store *A, x, b* data in memory visible to rank 0

- Library solves *Ax=b* using shared memory algorithms
- on the node

- App passes matrix and vector values to library data
- classes

- Using: SmpNode, SerialComm

## Example: Multicore + GPU Ax = b

- App
- Rank 0-3

- App
- Rank 4

- Lib
- Rank 0-3

- Lib
- Rank 4

- Kokkos objects allocate data according to concrete implementation.

- Library solves *Ax=b* using appropriate kernels

- App passes matrix and vector values to library data
- classes

- CPU
- memory

- GPU
- memory

- SMP
- kernel

- CUDA
- kernel

- Using: SmpNode, CudaNode, SerialComm, CudaComm

## Conclusion: Main Efforts

- Multi-precision support:
- Tpetra
- Tifpack
- Multi-precision Krylov in Belos.
- Multicore node support:
- Abstract node class.
- Compute kernels for multicore/GPUs.
- Concrete adapters for serial, Multicore, GPUs, ...
- Hybrid MPI-Only and MPI+shared memory runtime support.