Playa: Painless Linear Algebra

Version of the Day

Playa is a user-friendly, representation-independent system of components for development of high-performance parallel linear solvers, nonlinear solvers, and optimizers.

The five fundamental Playa object types are VectorType, VectorSpace, Vector, LinearOperator, and LinearSolver.

VectorType is an abstract factory that produces VectorSpace objects of a user-specified type, dimension, and distribution across processors.
VectorSpace is an abstract factory that produces Vector objects. This provides a consistent user interface for creating vectors without the client code needing to know anything about what sort of vector object is being created.
Vector objects represent vectors. In addition to the fundamental operations of vector addition and scalar-vector multiplication, arbitrary user-defined transformation and reduction operations can be implemented through application of templated functors. Some examples of such functors are in the PlayaFunctors namespace.
LinearOperator objects represent linear functions that map vector inputs to vector outputs.
LinearSolver is an object wrapper for algorithms for solving linear systems of equations $A{\bf x}={\bf b}.$ Solvers from several Trilinos packages are available through this interface: Amesos, AztecOO, and Belos, with preconditioning by Ifpack and ML. For specialized applications, user-defined solvers and preconditioners can be constructed out of Playa objects.

Contributors to Playa

Kevin Long -- Texas Tech University. Email kevin.long@ttu.edu
Paul Boggs -- Sandia National Laboratories, Livermore.
Victoria Howle -- Texas Tech University.
Robert Kirby -- Texas Tech University.

Quick start guide

Creating spaces and vectors

Vectors are rarely constructed directly by the user. Instead, they are usually created by a call to the createMember() method of VectorSpace. Similarly, a VectorSpace is usually constructed indirectly by one of the methods of a VectorType object. The user's choice of VectorType subclass specifies which internal object representation will be used.

EpetraVectorType specifies Epetra distributed vectors and matrices
SerialVectorType specifies serial vectors and dense matrices. If called in a SPMD program the spaces and vectors created will be replicated on each processor.

Here's how to create a Vector distributed with 10 elements per processer, represented internally by an Epetra_Vector object.

VectorType<double> vecType = new EpetraVectorType();
int numPerProc = 10;
VectorSpace<double> space = vecType.createEvenlyPartitionedSpace(MPIComm::world(), numPerProc);
Vector<double> x = space.createMember();

Spaces with more complicated distributions of data over processors can be created with the createSpace() method of VectorType.

Implicit operators

identityOperator() creates a matrix-free identity operator on a specified VectorSpace
zeroOperator() forms a matrix-free zero operator on specified range and domain spaces
diagonalOperator() forms a matrix-free diagonal operator with a specified vector
makeBlockOperator()
transposedOperator() creates an operator that performs a matrix-free transpose application
inverse() creates an implicit inverse of the specified LinearOperator with the solves carried out using the specified LinearSolver.

In addition, certain implicit operations can be specified using overloaded operators:

Operator addition and subtraction
```
LinearOperator<double> APlusB = A+B;
```
Scalar-operator multiplication
```
LinearOperator<double> aTimesA = a*A;
```
Operator composition
```
LinearOperator<double> AB = A*B;
```

All of these form implicit operators. For example, $z=ABx$

is computed by first computing $y=Bx$

, then

. The matrix-matrix product $ AB$

is never explicitly formed.

The implicit inverse and transpose operations can also be done through member functions of LinearOperator, for example,

LinearOperator<double> AInv = A.inverse(mySolver);
LinearOperator<double> ATrans = A.transpose();

Solvers

Linear solver algorithms are represented by subtypes of LinearSolver.

Preconditoners

The Preconditioner object stores a preconditioner with left and right operators. A preconditioner instance will rarely be created directly by the user; it will be created inside a solver by a user-specified subtype of PreconditionerFactory.

Solver Parameters

Linear solvers are usually created from Teuchos ParameterList objects.

LinearSolver<double> amesos = LinearSolverBuilder::createSolver("amesos.xml");

LinearSolver<double> solver = LinearSolverBuilder::createSolver(solverParams);

Amesos parameter example

Amesos is a serial sparse direct solver package.

<ParameterList>
  <ParameterList name="Linear Solver">
    <Parameter name="Type" type="string" value="Amesos"/>
    <Parameter name="Verbosity" type="int" value="0"/>
  </ParameterList>
</ParameterList>

Belos parameter examples

Parameters for Belos GMRES with ML preconditioning

<ParameterList>
  <ParameterList name="Linear Solver">
    <Parameter name="Type" type="string" value="Belos"/>
    <Parameter name="Maximum Iterations" type="int" value="1000"/>
    <Parameter name="Method" type="string" value="GMRES"/>
    <ParameterList name="Preconditioner">
       <Parameter name="Type" type="string" value="ML"/>
       <Parameter name="Problem Type" type="string" value="SA"/>
       <ParameterList name="ML Settings">
          <Parameter name="output" type="int" value="10"/>
       </ParameterList>
    </ParameterList> 
    <Parameter name="Convergence Tolerance" type="double" value="1e-12"/>
    <Parameter name="Output Frequency" type="int" value="0"/>
    <!-- num blocks is the restart size -->
    <Parameter name="Num Blocks" type="int" value="200"/>
    <Parameter name="Block Size" type="int" value="1"/>
  </ParameterList>
</ParameterList>

Parameters for Belos GMRES with ILU preconditioning

<ParameterList>
  <ParameterList name="Linear Solver">
    <Parameter name="Type" type="string" value="Belos"/>
    <Parameter name="Maximum Iterations" type="int" value="500"/>
    <Parameter name="Method" type="string" value="GMRES"/>
    <ParameterList name="Preconditioner">
       <Parameter name="Type" type="string" value="Ifpack"/>
       <Parameter name="Prec Type" type="string" value="ILU"/>
       <Parameter name="Overlap" type="int" value="1"/>
       <ParameterList name="Ifpack Settings">
          <Parameter name="fact: level-of-fill" type="int" value="2"/>
       </ParameterList>
    </ParameterList> 
    <Parameter name="Convergence Tolerance" type="double" value="1e-10"/>
    <Parameter name="Output Frequency" type="int" value="0"/>
    <!-- num blocks is the restart size -->
    <Parameter name="Num Blocks" type="int" value="200"/>
    <Parameter name="Block Size" type="int" value="1"/>
    <!-- set verbosity to 33 or 41 for useful diagnostics -->
    <Parameter name="Verbosity" type="int" value="0"/>
  </ParameterList>
</ParameterList>

Aztec parameter examples

<ParameterList>
  <ParameterList name="Linear Solver">
    <Parameter name="Max Iterations" type="int" value="500"/>
    <Parameter name="Method" type="string" value="GMRES"/>
    <ParameterList name="Preconditioner">
       <Parameter name="Type" type="string" value="ML"/>
       <Parameter name="Problem Type" type="string" value="SA"/>
       <ParameterList name="ML Settings">
          <Parameter name="output" type="int" value="10"/>
       </ParameterList>
    </ParameterList> 
    <Parameter name="Tolerance" type="double" value="1e-12"/>
    <Parameter name="Type" type="string" value="Aztec"/>
    <Parameter name="Verbosity" type="int" value="0"/>
  </ParameterList>
</ParameterList>