# Presentation: 850340

## An Overview of TSFCore

**Roscoe A. Bartlett**
** **
**9211, Optimization and Uncertainty Estimation**

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

## Get most recent copy at:  Trilinos/doc/TSFCore

## Nonlinear Equations : Foundation for all our Work!

- _Applications_
-   Discretized PDEs (e.g. finite element, finite volume, finite difference etc.)
-   Network problems (e.g. Xyce)

## Nonlinear Equations : Sensitivities

**Related Algorithms**
-   Gradient-based optimization
-   SAND
-   NAND
-   Nonlinear equations (NLS)
-   Multidisciplinary analysis
-   Linear (matrix) analysis
-   Block iterative solvers
-   Eigenvalue problems
-   Uncertainty quantification
-   SFE
-   Stability analysis / continuation
-   Transients (ODEs, DAEs)

- B. van Bloemen Waanders, R. A. Bartlett, K. R. Long and P. T. Boggs. Large Scale Non-Linear Programming: PDE Applications and Dynamical Systems, Sandia National Laboratories, SAND2002-3198, 2002

## Applications, Algorithms, Linear-Algebra Software

**APP : Application (e.g. MPSalsa, Xyce, SIERRA, NEVADA etc.)**
**LAL : Linear-Algebra Library (e.g. Petra/Ifpack, PETSc, Aztec etc.)**
**ANA : Abstract Numerical Algorithm (e.g. optimization, nonlinear solvers, stability analysis, SFE, transient solvers etc.)**

**Key points**
-   Complex algorithms
-   Complex software
-   Complex interfaces
-   Complex computers
-   Duplication of effort?

## TSFCore

**Key points**
-   Maximizing development impact
-   Software can be run on more sophisticated computers
-   Fosters improved algorithm development

## Requirements for TSFCore

- TSFCore should:
- Be portable to ASCI platforms
- Provide for stable and accurate numerical computations
- Represent a minimal but complete interface that will result in implementations that are:
- Near optimal in computational speed
- Near optimal in storage
- Be independent of computing environment (SPMD, MS, CS etc.)
- Be easy to develop adapters for existing libraries (e.g. Epetra, PETSc etc.)

## Example ANA : Linear Conjugate Gradient Solver

## TSFCore : Basic Linear Algebra Interfaces

**Warning! Unified Modeling Langage (UML) Notation!**

- <<create>>

## TSFCore : Basic Linear Algebra Interfaces

- <<create>>

## TSFCore : Basic Linear Algebra Interfaces

- <<create>>

## TSFCore : Basic Linear Algebra Interfaces

- <<create>>

## TSFCore : Basic Linear Algebra Interfaces

**The Key to success!**
**Reduction/Transformation Operators**** **
-   Supports all needed vector operations
-   Data/parallel independence
-   Optimal performance

- R. A. Bartlett, B. G. van Bloemen Waanders and M. A. Heroux. Vector Reduction/Transformation Operators, Accepted to ACM TOMS, 2003

## Background for TSFCore

- 1996 : Hilbert Class Library (HCL), [Symes and Gockenbach]
- Abstract vector spaces, vectors, linear operators
- 2000 : Epetra, [Heroux]
- Concrete multi-vectors
- 2001 : Trilinos Solver Framework (TSF) 0.1, [Long]
- 2001 : AbstractLinAlgPack (ALAP) (MOOCHO LA interfaces), [Bartlett]
- Reduction/transformation operators (RTOp)
- Abstract multi-vectors

## TSFCore : Basic Linear Algebra Interfaces

## TSFCore Details

- All interfaces are templated on **Scalar** type (support real and complex)
- Smart reference counted pointer class **Teuchos::RefCountPtr<>** used for all dynamic memory management
- Many operations have default implementations based on very few pure virtual methods
**RTOp** operators (and wrapper functions) are provided for many common level-1 vector and multi-vector operations
- Default implementation provided for **MultiVector** (**MultiVectorCols**)
- Default implementations provided for serial computation: **VectorSpace** (**SerialVectorSpace**), **VectorSpaceFactory** (**SerialVectorSpaceFactory**), **Vector** (**SerialVector**)

## Vector-Vector Operations Provided with TSFCore

- namespace TSFCore {
-     template<class Scalar> Scalar sum( const Vector<Scalar>& v );                   // result = sum(v(i))
-     template<class Scalar> Scalar norm_1( const Vector<Scalar>& v );                // result = ||v||1
-     template<class Scalar> Scalar norm_2( const Vector<Scalar>& v );                // result = ||v||2
-     template<class Scalar> Scalar norm_inf( const Vector<Scalar>& v_rhs );          // result = ||v||inf
-     template<class Scalar> Scalar dot( const Vector<Scalar>& x
-                                        ,const Vector<Scalar>& y );                  // result = x'*y
-     template<class Scalar> Scalar get_ele( const Vector<Scalar>& v, Index i );      // result = v(i)
-     template<class Scalar> void set_ele( Index i, Scalar alpha
-                                          ,Vector<Scalar>* v );                      // v(i) = alpha
-     template<class Scalar> void assign( Vector<Scalar>* y, const Scalar& alpha );   // y = alpha
-     template<class Scalar> void assign( Vector<Scalar>* y
-                                         ,const Vector<Scalar>& x );                 // y = x
-     template<class Scalar> void Vp_S( Vector<Scalar>* y, const Scalar& alpha );     // y += alpha
-     template<class Scalar> void Vt_S( Vector<Scalar>* y, const Scalar& alpha );     // y *= alpha
-     template<class Scalar> void Vp_StV( Vector<Scalar>* y, const Scalar& alpha
-                                         ,const Vector<Scalar>& x );                 // y = alpha*x + y
-     template<class Scalar> void ele_wise_prod( const Scalar& alpha
-         ,const Vector<Scalar>& x, const Vector<Scalar>& v, Vector<Scalar>* y );     // y(i)+=alpha*x(i)*v(i)
-     template<class Scalar> void ele_wise_divide( const Scalar& alpha
-         ,const Vector<Scalar>& x, const Vector<Scalar>& v, Vector<Scalar>* y );     // y(i)=alpha*x(i)/v(i)
-     template<class Scalar> void seed_randomize( unsigned int );                     // Seed for randomize()
-     template<class Scalar> void randomize( Scalar l, Scalar u, Vector<Scalar>* v ); // v(i) = random(l,u)
- } // end namespace TSFCore

## TSFCore : Vectors and Vector Spaces

**TSFCore : Vectors and Vector Spaces**

**C++ code:**
- template<class Scalar>
- Scalar foo( const VectorSpace<Scalar>& S )
- {
-     Teuchos::RefCountPtr<Vector<Scalar> >
-         x = S.createMember(),                      // create x
-         y = S.createMember();                      // create y
-     assign( &*x, 1.0 );                            // x = 1
-     randomize( -1.0, +1.0, &*y );                  // y = rand(-1,1)
-     Vp_StV( &*y, -2.0, *x );                       // y += -2.0 * x
-     Scalar gamma = dot(*x,*y);                     // gamma = x’*y
-     return gamma;
- }

**Mathematical notation:**

## TSFCore : Applying a Linear Operator

**TSFCore : Applying a Linear Operator**

_**C++ Prototype:**_
- namespace TSFCore {
-     enum ETransp { NOTRANS, TRANS, CONJTRANS };
-     template<class Scalar>
-     class LinearOp : public virtual OpBase<Scalar> {
-     public:
-     virtual void apply(
-             ETransp M_trans, const Vector<Scalar> &x, Vector<Scalar> *y
-             ,Scalar alpha = 1.0, Scalar beta = 0.0
-             ) const = 0;
-     };
- }

_**Example:**_
- template<class Scalar>
- void myOp( const Vector<Scalar> &x, const LinearOp<Scalar> &M
-            ,Vector<Scalar> *y )
- {
-     M.apply( NOTRANS, x, y );
- }

## Example ANA : Linear Conjugate Gradient Solver

## Multi-vector Conjugate-Gradient Solver : Single Iteration

**Multi-vector Conjugate-Gradient Solver : Single Iteration**

- template<class Scalar>
- void CGSolver<Scalar>::doIteration(
-     const LinearOp<Scalar> &M, ETransp opM_notrans
-     ,ETransp opM_trans, MultiVector<Scalar> *X, Scalar a
-     ,const LinearOp<Scalar> *M_tilde_inv
-     ,ETransp opM_tilde_inv_notrans, ETransp opM_tilde_inv_trans
-     ) const
- {
-     const Index m = currNumSystems_;
-     int j;
-     if( M_tilde_inv )
-         M_tilde_inv->apply( opM_tilde_inv_notrans, *R_, &*Z_ );
-     else
-         assign( &*Z_, *R_ );
-     dot( *Z_, *R_, &rho_[0] );
-     if( currIteration_ == 1 ) {
-         assign( &*P_, *Z_ );
-     }
-     else {
-         for(j=0;j<m;++j) beta_[j] = rho_[j]/rho_old_[j];
-         update( *Z_, &beta_[0], 1.0, &*P_ );
-     }
-     M.apply( opM_notrans, *P_, &*Q_ );
-     dot( *P_, *Q_, &gamma_[0] );
-     for(j=0;j<m;++j) alpha_[j] = rho_[j]/gamma_[j];
-     update( &alpha_[0], +1.0, *P_, X )
-     update( &alpha_[0], -1.0, *Q_, &*R_ );
- }

## The TSFCore Trilinos package

- packages/TSFCore
- src
- interfaces
- Core : VectorSpace, Vector, LinearOp etc ...
- Solvers : Iterative linear solver interfaces (unofficial!)
- Nonlin : Nonlinear problem interfaces (unofficial!)
- utilities
- Core : Testing etc  ...
- Solvers : Some iterative solvers (CG, BiCG, GMRES)
- Nonlin : Testing etc ...
- adapters
- mpi-base : Node classes for MPI-based vector spaces
- Epetra : EpetraVectorSpace, EpetraVector etc ...
- examples ...

## TSFCore::Nonlin : Interfaces to Nonlinear Problems

**Supported Areas**
-   NAND optimization
-   SAND optimization
-   Nonlinear equations
-   Multidisciplinary analysis
-   Stability analysis / continuation
-   SFE

**Notes:**
Bookman old style is the closet to the latex font.

## TSFCore::Nonlin : Interfaces to Nonlinear Problems

**response functions**

**TSFCore::Nonlin : Interfaces to Nonlinear Problems**

**Supported Areas**
-   SAND
-   Nonlinear equations
-   Multidisciplinary analysis
-   Stability analysis / continuation
-   SFE

## Summary

_**SAND Reports**_
- R. A. Bartlett, M. A. Heroux and K. R. Long.  TSFCore : A Package of Light-Weight Object-Oriented Abstractions for the Development of Abstract Numerical Algorithms and Interfacing to Linear Algebra Libraries and Applications, Sandia National Laboratories, SAND2003-1378, 2003
- R. A. Bartlett, TSFCore::Nonlin : An Extension of TSFCore for the Development of Nonlinear Abstract Numerical Algorithms and Interfacing to Nonlinear Applications, Sandia National Laboratories, SAND2003-1377, 2003
**Location:** Trilinos/doc/TSFCore

**Notes:**
Bookman old style is the closet to the latex font.

## The End

**Thank You!**

## Extra Slides

## Examples of Non-Standard Vector Operations

- Currently in MOOCHO :
-    > 40 vector operations!

## Goals for a Vector Interface

- Compute efficiency                                            =>  Near optimal performance
- Optimization developers add new operations    => Independence of linear algebra        .                                                                         library developers
- Compute environment independence                 => Flexible optimization software
- Minimal number of methods                              => Easy to write adapters

## (1)  Linear algebra library allows direct access to vector elements

- (2)  Optimizer-specific interfaces
- (3) General-purpose primitive vector operations

##   Communication between sets of subvectors only for  NULL, oprr( 1 , 2 )  2

**Reduction/Transformation Operators (RTOp) Defined**
*z** *1*i *... *z** q**i*   opt( *i* , *v** *1*i *... *v** p**i *,* z** *1*i *... *z** q**i *)        element-wise transformation
-   opr( *i *, *v** *1*i *... *v** p**i *,* z** *1*i *... *z** q**i* )        element-wise reduction
- 2              oprr( 1 , 2 )                                 reduction of intermediate reduction objects
*  **v** *1* *... *v** p* ***R*** n  : *p* non-mutable input vectors
*  **z** *1* *... *z** q** ****R*** n   : *q* mutable input/output vectors
-   : reduction target object (many be  non-scalar (e.g. {*y**k *,*k*}), or NULL)

**Key to Optimal Performance**
-   opt(...) and opr(...) applied to entire sets of subvectors (*i *= *a*...*b*) independently:
*z** *1*a:b *... *z** q**a:b **, *   op( *a*,* b* , *v** *1*a:b *... *v** p**a:b *,* z** *1*a:b *... *z** q**a:b** , ** *)
-   Communication between sets of subvectors only for  NULL, oprr( 1 , 2 )  2

## Object-Oriented Design for User Defined RTOp Operators

**Advantages:**
-   Functionality
-   Linear-algebra implementations can be changed with no impact on optimizer
-   Optimizer developers can unilaterally add new vector operations
-   Performance
-   Near optimal performance (large subvectors)
-   Multiple simultaneous global reductions => no sequential bottlenecks
-   No unnecessary temporary vectors or multiple vector read/writes
-   **Disadvantages:**
-   New concepts, initially harder to understand interfaces?

## RTOp vs. Primitives : Communication

**Compare**
**RTOp (all-at-once reduction (i.e. ISIS++ QMR solver))**
- {, , , ,  }  { (*x*T *x*)1/2, (*v*T *v*)1/2, (*w*T *w*)1/2, *w*T *v*, *v*T *t* }
**Primitives (5 separate reductions)**
-  (*x*T *x*)1/2,  (*v*T *v*)1/2,  (*w*T *w*)1/2,  *w*T *v*,  *v*T

- * 128 processors on CPlant®

## RTOp vs. Primitives : Multiple Ops and Temporaries

**Compare**
**RTOp (all-at-once reduction)**
- { max ** : *x + **** d *** } = min{ max( (** - *x**i*)/*d**i*, *0* ), for *i = 1 ... n* }  **
**Primitives (5 temporaries, 6 vector operations)**
- -*x**i*  *u**i*,   *x**i* + **  *v**i*,   *v**i** */ *d**i*  *w**i*,   *0*  *y**i*,   max{*w**i*,*y**i*}  *z**i*,   min{*z**i*,*i=1...n*}  

- * 1 processor (gcc 3.1 under Linux)

## Question:  Does OO C++ allow for good scalability for massively parallel computing (i.e. 100 to 10000 processors)?

**Parallel Scalability of MOOCHO**

**Where is the parallel bottleneck?**
**Is it OO C++ or MPI?**

*** Red Hat Linux cluster (4 nodes)**
** **2.0 GHz Intel P4 processors
-  MPICH 1.2.2.1

**                                        ****Answer => MPI**

**Serial overhead of MOOCHO (n=2, N****p****=1)**
**       **0.41 milliseconds per rSQP iteration

**Overhead of MPI communication (N****p****=4)**
**       ** 0.42  milliseconds per global reduction