# Presentation: 009438

## Trilinos 101 (Part II)Creating and managing linear algebra data in Trilinos: Data management using Epetra and Teuchos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:**
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

- Creating Objects.
- 2D Objects.
- Teuchos tidbits.
- Performance Optimizations.

## A Simple Epetra/AztecOO Program

- int main(int argc, char *argv[]) {
-   MPI_Init(&argc,&argv); // Initialize MPI, MpiComm
-   Epetra_MpiComm Comm( MPI_COMM_WORLD );

**A Simple Epetra/AztecOO Program**

-   // ***** Create x and b vectors *****
-   Epetra_Vector x(Map);
-   Epetra_Vector b(Map);
-   b.Random(); // Fill RHS with random #s

- // ***** Create an Epetra_Matrix  tridiag(-1,2,-1) *****
-   Epetra_CrsMatrix A(Copy, Map, 3);
-   double negOne = -1.0; double posTwo = 2.0;
-   for (int i=0; i<NumMyElements; i++) {
-     int GlobalRow = A.GRID(i);
-     int RowLess1 = GlobalRow - 1;
-     int RowPlus1 = GlobalRow + 1;
-     if (RowLess1!=-1)
-        A.InsertGlobalValues(GlobalRow, 1, &negOne, &RowLess1);
-     if (RowPlus1!=NumGlobalElements)
-        A.InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus1);
-     A.InsertGlobalValues(GlobalRow, 1, &posTwo, &GlobalRow);
-   }
- A.FillComplete(); // Transform from GIDs to LIDs

- // ***** Map puts same number of equations on each pe *****
-   int NumMyElements = 1000 ;
-   Epetra_Map Map(-1, NumMyElements, 0, Comm);
-   int NumGlobalElements = Map.NumGlobalElements();

- // ***** Report results, finish ***********************
-   cout << "Solver performed " << solver.NumIters()
-           << " iterations." << endl
-           << "Norm of true residual = "
-           << solver.TrueResidual()
-           << endl;
-   MPI_Finalize() ;
-   return 0;
- }

-   // ***** Create/define AztecOO instance, solve *****
-   AztecOO solver(problem);
-   solver.SetAztecOption(AZ_precond, AZ_Jacobi);
-   solver.Iterate(1000, 1.0E-8);

- // ***** Create Linear Problem *****
-   Epetra_LinearProblem problem(&A, &x, &b);

- // Header files omitted...
- int main(int argc, char *argv[]) {
- Epetra_SerialComm Comm();

## Typical Flow of Epetra Object Construction

- Construct Comm

- Construct Map

- Construct *x*

- Construct *b*

- Construct *A*

-  Any number of Comm objects can exist.
-  Comms can be nested (e.g., serial within MPI).

-  Maps describe parallel layout.
-  Maps typically associated with more than one comp object.
-  Two maps (source and target) define an export/import object.

-  Computational objects.
-  Compatibility assured via common map.

## Trilinos/packages/epetraext/test/inout/ (Part 1)

- EpetraExt::BlockMapToMatrixMarketFile("Test_map.mm", *map, "Official EpetraExt test map",
-                    "This is the official EpetraExt test map generated by the EpetraExt regression tests");
- EpetraExt::RowMatrixToMatrixMarketFile("Test_A.mm", *A, "Official EpetraExt test matrix",
-                    "This is the official EpetraExt test matrix generated by the EpetraExt regression tests");
- EpetraExt::VectorToMatrixMarketFile("Test_x.mm", *x, "Official EpetraExt test initial guess",
-                    "This is the official EpetraExt test initial guess generated by the EpetraExt regression tests");
- EpetraExt::VectorToMatrixMarketFile("Test_xexact.mm", *xexact, "Official EpetraExt test exact solution",
-                    "This is the official EpetraExt test exact solution generated by the EpetraExt regression tests");
- EpetraExt::VectorToMatrixMarketFile("Test_b.mm", *b, "Official EpetraExt test right hand side",
-                    "This is the official EpetraExt test right hand side generated by the EpetraExt regression tests");

## Trilinos/packages/epetraext/test/inout/ (part 2)

- Epetra_Map * map1;
- Epetra_CrsMatrix * A1;
- Epetra_Vector * x1;
- Epetra_Vector * b1;
- Epetra_Vector * xexact1;
- EpetraExt::MatrixMarketFileToMap("Test_map.mm", comm, map1);
- if (map->SameAs(*map1))
-   if (verbose) cout << "Maps are equal.  In/Out works." << endl;
- else
-   if (verbose) cout << "Maps are not equal.  In/Out fails." << endl;
- // Read Matrix:
- EpetraExt::MatrixMarketFileToCrsMatrix("Test_A.mm", *map1, A1);
- // Alternates:
- EpetraExt::MatrixMarketFileToCrsMatrix("Test_A.mm", comm, A1); // Only need comm
- EpetraExt::MatlabFileToCrsMatrix("Test_A.dat", comm, A1); // Read (row,col,val) (no MM header).
- // Read Vector (Equivalent MultiVector versions.
- EpetraExt::MatrixMarketFileToVector("Test_x.mm", *map1, x1);
- EpetraExt::MatrixMarketFileToVector("Test_xexact.mm", *map1, xexact1);
- EpetraExt::MatrixMarketFileToVector("Test_b.mm", *map1, b1);

## Inout Summary

- Reads from/writes to Matlab or Matrix Market compatible files.
- Works for any distributed map, matrix, vector or multivector.
- If map is read in on the same number of processor it was written to, the map layout will be preserved.

## Details about Epetra Maps

- Note:  Focus on Maps (not BlockMaps).
- Getting beyond standard use case...

## 1-to-1 Maps

*1-to-1 map* (defn):  A map is 1-to-1 if each GID appears only once in the map (and is therefore associated with only a single processor).
- Certain operations in parallel data repartitioning require 1-to-1 maps.  Specifically:
- The source map of an import must be 1-to-1.
- The target map of an export must be 1-to-1.
- The domain map of a 2D object must be 1-to-1.
-  The range map of a 2D object must be 1-to-1.

## 2D Objects: Four Maps

- Epetra 2D objects:
- CrsMatrix
- CrsGraph
- VbrMatrix
- Have four maps:
**RowMap**: On each processor, the GIDs of the **rows** that processor will “manage”.
**ColMap**: On each processor, the GIDs of the **columns** that processor will “manage”.
**DomainMap**: The layout of domain objects 	(the *x* vector/multivector in *y=Ax*).
**RangeMap**: The layout of range objects 	(the *y* vector/multivector in *y=Ax*).

- Must be 1-to-1 maps!!!

- Typically a 1-to-1 map

- Typically **NOT **a 1-to-1 map

## Sample Problem

- =

*y*

*A*

*x*

## Case 1: Standard Approach

- RowMap 	= {0, 1}
- ColMap 	= {0, 1, 2}
- DomainMap 	= {0, 1}
- RangeMap 	= {0, 1}

- First 2 rows of *A*, elements of *y* and elements of *x*, kept on PE 0.
- Last row of *A*, element of *y* and element of *x*, kept on PE 1.

**PE 0 Contents**

**PE 1 Contents**

- RowMap 	= {2}
- ColMap 	= {1, 2}
- DomainMap 	= {2}
- RangeMap 	= {2}

- Notes:
- Rows are wholly owned.
- RowMap=DomainMap=RangeMap (all 1-to-1).
- ColMap is NOT 1-to-1.
- Call to FillComplete: **A.FillComplete(); // Assumes**

- =

*y*

*A*

*x*

- Original Problem

## Case 2: Twist 1

**Case 2: Twist 1**

- RowMap 	= {0, 1}
- ColMap 	= {0, 1, 2}
- DomainMap 	= {1, 2}
- RangeMap 	= {0}

- First 2 rows of *A*, first element of *y* and last 2 elements of *x*, kept on PE 0.
- Last row of *A*, last 2 element of *y* and first element of *x*, kept on PE 1.

**PE 0 Contents**

**PE 1 Contents**

- RowMap 	= {2}
- ColMap 	= {1, 2}
- DomainMap 	= {0}
- RangeMap 	= {1, 2}

- Notes:
- Rows are wholly owned.
- RowMap is NOT = DomainMap 	     is NOT = RangeMap (all 1-to-1).
- ColMap is NOT 1-to-1.
- Call to FillComplete: **A.FillComplete(DomainMap, RangeMap);**

- =

*y*

*A*

*x*

- Original Problem

## Case 2: Twist 2

- RowMap 	= {0, 1}
- ColMap 	= {0, 1}
- DomainMap 	= {1, 2}
- RangeMap 	= {0}

- First row of *A*, part of second row of *A, *first element of *y* and last 2 elements of *x*, kept on PE 0.
- Last row, part of second row of *A, *last 2 element of *y* and first element of *x*, kept on PE 1.

**PE 0 Contents**

**PE 1 Contents**

- RowMap 	= {1, 2}
- ColMap 	= {1, 2}
- DomainMap 	= {0}
- RangeMap 	= {1, 2}

- Notes:
- Rows are NOT wholly owned.
- RowMap is NOT = DomainMap 	     is NOT = RangeMap (all 1-to-1).
- RowMap and ColMap are NOT 1-to-1.
- Call to FillComplete: **A.FillComplete(DomainMap, RangeMap);**

- =

*y*

*A*

*x*

- Original Problem

## What does FillComplete Do?

- A bunch of stuff.
- One task is to create (if needed) import/export objects to support distributed matrix-vector multiplication:
- If ColMap  DomainMap, create Import object.
- If RowMap  RangeMap, create Export object.
- A few rules:
- Rectangular matrices will *always* require:
- A.FillComplete(DomainMap,RangeMap);
- DomainMap and RangeMap *must be 1-to-1.*

## Teuchos::RefCountPtrGoal: Empty Destructors

- Goal of memory allocation: Empty destructors.
- RefCountPtr important for this goal.
- Sample attribute list:

- private:
-   Epetra_Map cmsMap_;
-   Epetra_Map densityMap_;
-   Epetra_Map block2Map_;
-   Teuchos::RefCountPtr<Epetra_CrsMatrix> cmsOnDensityMatrix_;
-   Teuchos::RefCountPtr<Epetra_Vector> cmsOnCmsMatrix_;
-   Teuchos::RefCountPtr<Epetra_Vector> densityOnDensityMatrix_;
-   Teuchos::RefCountPtr<Epetra_Vector> densityOnCmsMatrix_;
-   Epetra_IntSerialDenseVector indices_;
-   Epetra_SerialDenseVector values_;
-   Ifpack factory_;
- };
- #endif /* DFT_POLYA22_EPETRA_OPERATOR_H */

## Exception Testing

- Teuchos exception macros: 	Standard approach to handling errors.

- int dft_PolyA22_Epetra_Operator::ApplyInverse(const Epetra_MultiVector& X,
- 				                        Epetra_MultiVector& Y) const {
-   TEST_FOR_EXCEPT(!X.Map().SameAs(OperatorDomainMap()));
-   TEST_FOR_EXCEPT(!Y.Map().SameAs(OperatorRangeMap()));
-   TEST_FOR_EXCEPT(Y.NumVectors()!=X.NumVectors());

## Epetra Performance Optimization GuideSAND2005-1668

- Topics:
- 3rd Party Libraries: BLAS and LAPACK
- Epetra_MultiVector Data Layout
- Epetra_CrsGraph Construction
- Epetra_CrsMatrix Construction
- Selecting the Right Sparse Matrix Class
- Parallel Data Redistribution
- General Practices
- Tiered Approach to practices.

## Practice Categories

**VSR**

**SR**

**R**

- Each practice falls into one of three categories:

**Very Strongly Recommended** - Practices necessary for Epetra to perform well.

**Strongly Recommended **** ** - Practices that are definitely a good thing or that have proved to be valuable.

**Recommended **- Practices that are probably a good idea.

## Epetra Computational Classes

- Epetra contains base class Epetra_CompObject.
- Small class, but classes that derive from it are exactly those which are the focus of guide.

## Fast BLAS

- Link to high-performance BLAS for Epetra_MultiVector, Epetra_VbrMatrix and Epetra_SerialDense performance
- Options:
- Atlas (All platforms).
- GOTO BLAS (PC).
- vecLib (Mac).
- Native BLAS on legacy platforms.

**VSR**

## Epetra_MultiVector Layout

- Create Epetra_MultiVector objects with strided storage.
- Improved block vector updates, block dot products.
- Improve sparse matrix multiplication/solves.

**VSR**

## Epetra_CrsGraph/Matrix Construction 

**Construct Epetra_CrsGraph objects first**:
- When constructing multiple Epetra_CrsMatrix or Epetra_VbrMatrix objects, much of the overhead in sparse matrix construction is related solely to the graph structure of the matrix.
- By pre-constructing the Epetra_CrsGraph object, this expense is incurred only once and then amortized over multiple matrix constructions.
- Note: Even when constructing a single matrix, it is often the case that matrix construction is done within a nonlinear iteration or a time-stepping iteration.  In both of these cases, it is best to construct the Epetra_CrsGraph object one time and then reuse it.
**When constructing Epetra_CrsGraph/matrix objects, carefully manage the nonzero profile and set StaticProfile to ‘true’ :**
- Although it is very convenient to use the flexible insertion capabilities of Epetra_CrsGraph/Matrix, the performance and memory cost can be substantial.
- StaticProfile is an optional argument to the Epetra_CrsGraph/Matrix constructors that forces the constructor to pre-allocate all storage for the graph, using the argument NumIndicesPerRow as a strict upper limit on the number of indices that will be inserted.
**After calling FillComplete(), call OptimizeStorage() : **
- The OptimizeStorage() method frees memory that is used only when submitting data to the object.
- Also, the storage that remains is packed for better cache memory performance due to better spatial locality.
- If StaticProfile is true, the packing of data is fast and cheap, if needed at all .
- If StaticProfile is false, then the graph data is not contiguously stored, so we must allocate contiguous space for all entries, copy data from the existing arrays, and then delete the old arrays.  This can be a substantial overhead in memory costs and could even double the high-water memory mark.  However, the performance improvement for matrix operations can be 20% to a full doubling of performance.

**VSR**

**VSR**

**VSR**

## Epetra_CrsGraph/Matrix (cont)

**Use Epetra_FECrsMatrix if you have many repeated indices.  Then extract the associated Epetra_CrsGraph for subsequent use: **
- Although there is no Epetra_FECrsGraph class (something that we may introduce in the future), it is typically best to use the Epetra_FECrsMatrix class to construct matrices when you have many repeated indices.
- This is typically the case when forming global stiffness matrices from local element stiffness matrices in a finite element assembly loop.
- Note:  Even though there is no Epetra_FECrsGraph class, once an Epetra_FECrsMatrix has been constructed, call it myFEMatrix, there is an Epetra_CrsGraph that can be accessed via myFEMatrix.Graph()  and used to construct further matrices that have the same pattern as myFEMatrix.
**When inserting indices into Epetra_CrsGraph/Matrix objects, avoid large numbers of repeated indices:**
- To reduce time costs, indices that are inserted into a row are not checked for redundancy at the time they are inserted.  Instead, when FillComplete() is called, indices are sorted and then redundant indices are removed.
- Because of this approach, repeated indices increase the memory used by a given row, at least until FillComplete() is called.
**Submit large numbers of graph/matrix indices at once: **
- When inserting indices, submit many entries at once if possible.
- Ideally, it is best to submit all indices for a given row at once.

**VSR**

**SR**

**SR**

## Picking a Sparse Matrix Class

**Epetra_CrsMatrix**: Compressed Row Storage scalar entry matrix.
**Epetra_FECrsMatrix**: Finite Element Compressed Row Storage matrix.
**Epetra_VbrMatrix**: Variable Block Row block entry matrix.
**Epetra_FEVbrMatrix**: Finite Element Variable Block Row matrix.
- Note: **Epetra_JadOperator** does not implement Epetra_RowMatrix, but is worth  mentioning for users of vector processors.  This class is a first version to support vectorization of sparse matrix-vector multiplication.

## Picking a Sparse Matrix Class (cont)

**Use Epetra_FEVbrMatrix to construct a matrix with dense block entries and repeated summations into block entries**
**Use Epetra_VbrMatrix to construct a matrix with dense block entries and few repeated submissions:**
- Block entries should be of size 4 or larger before Vbr formats are preferable to Crs formats.
**Use Epetra_FECrsMatrix to construct a matrix with scalar entries and repeated summations into entries: **
- Epetra_FECrsMatrix is designed to handle problems such as finite element, finite volume or finite difference problems where a single degree of freedom is being tracked at a single mesh point, and the matrix entries are being summed into the global stiffness matrix one finite element or control volume cell at a time.
**Use Epetra_CrsMatrix in all other cases:**
- Epetra_CrsMatrix is the simplest and most natural matrix data structure for people who are used to thinking about sparse matrices.  Unless you have one of the above three situations, you should use Epetra_CrsMatrix.
**Use Epetra_RowMatrix to access matrix entries: **
- If you are writing code to use an existing Epetra sparse matrix object, use the Epetra_RowMatrix interface to access the matrix.
- This will provide you compatibility with all of Epetra’s sparse matrix classes, and allow users of your code to provide their own custom implementation of Epetra_RowMatrix if needed.

**VSR**

**R**

**VSR**

**VSR**

**VSR**

## Other RowMatrix adapters

- In addition, there are numerous other implementation of Epetra_RowMatrix provided in other packages.  For example:
- Epetra_MsrMatrix:
- This class is provided within AztecOO.
- The constructor for Epetra_MsrMatrix takes a single argument, namely an existing AZ_DMSR matrix.
- Epetra_MsrMatrix does not make a copy of the Aztec matrix.  Instead it make the Aztec matrix act like an Epetra_RowMatrix object.
- Ifpack Filters:
- Ifpack uses the Epetra_RowMatrix interface to provide modified views of existing Epetra_RowMatrix objects.
-  For example, Ifpack_DropFilter creates a new Epetra_RowMatrix from an existing one by dropping all matrix values below a certain tolerance.

## Parallel Data Redistribution

- Use EpetraExt to balance load:
- Although an advanced feature, any serious use of Epetra for scalable performance must ensure that work and data are balanced across the parallel machine.
- EpetraExt provides an interface to the Zoltan library that can greatly improve the load balance for typical sparse matrix computations.
- Use Epetra_OffsetIndex to improve performance of multiple redistributions:
- Often the data distribution that is optimal for constructing sparse matrices is not optimal for solving the related system of equations.
- As a result, the matrix will be constructed in one distribution, and then redistributed for the solve.
- Epetra_OffsetIndex precomputes the offsets for the redistribution pattern, making repeated redistributions cheaper.

## General Practices

- Check method return codes:
- Almost all methods of Epetra classes provide an integer return argument.
- This argument is set to zero if no errors or warning occurred.
- However, a number of performance-sensitive methods will communicate potential performance problems by returning a positive error code.
- For example, if the matrix insertion routine must re-allocate memory because the allocated row storage is too small, a return code of 1 is set.  Since memory reallocation can be expensive, this information can be an indication that the user should increase the nonzero profile values (NumEntriesPerRow).
- Compile Epetra with aggressive optimization:
- Epetra performance can benefit greatly from aggressive compiler optimization settings.
- In particular, aggressive optimization for FORTRAN can be especially important. Configuring Epetra (or Epetra as part of Trilinos) with the following compiler flags works well on 32-bit PC platforms with the GCC compilers:
- ../configure CXXFLAGS="-O3" CFLAGS="-O3" \
- FFLAGS="-O5 -funroll-all-loops -malign-double"

## Summary

- 2D objects:
- Epetra provides very flexible support for:
- Placing any matrix entry on any processor, and
- Even sharing a matrix entry across processors.
- Performance Optimizations Matter.