slide1 l.
Download
Skip this Video
Loading SlideShow in 5 Seconds..
Next Generation Iterative Solvers: Belos & Anasazi November 3 rd , 9:30-10:30am Mike Heroux Ulrich Hetmaniuk Rich L PowerPoint Presentation
Download Presentation
Next Generation Iterative Solvers: Belos & Anasazi November 3 rd , 9:30-10:30am Mike Heroux Ulrich Hetmaniuk Rich L

Loading in 2 Seconds...

play fullscreen
1 / 25

Next Generation Iterative Solvers: Belos & Anasazi November 3 rd , 9:30-10:30am Mike Heroux Ulrich Hetmaniuk Rich L - PowerPoint PPT Presentation


  • 158 Views
  • Uploaded on

Next Generation Iterative Solvers: Belos & Anasazi November 3 rd , 9:30-10:30am Mike Heroux Ulrich Hetmaniuk Rich Lehoucq Heidi Thornquist.

loader
I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.
capcha
Download Presentation

PowerPoint Slideshow about 'Next Generation Iterative Solvers: Belos & Anasazi November 3 rd , 9:30-10:30am Mike Heroux Ulrich Hetmaniuk Rich L' - ranger


An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.


- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript
slide1

Next Generation Iterative Solvers: Belos & AnasaziNovember 3rd, 9:30-10:30amMike HerouxUlrich HetmaniukRich LehoucqHeidi Thornquist

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.

outline
Outline
  • Background
    • AztecOO and ARPACK
    • SLEPc
  • Iterative solver design requirements
    • Generic programming
  • Belos and Anasazi
    • What’s in Trilinos 4.0
    • Current Interfaces
    • Future Development
  • Concluding remarks
background
Background
  • Several iterative linear solver / eigensolver libraries exist:
    • PETSc, SLAP, LINSOL, Aztec(OO), …
    • SLEPc, LOBPCG (hypre), ARPACK, …
  • None are (completely) written in C++
  • None of the linear solver libraries can efficiently deal with multiple right-hand sides
  • Stopping criteria are predetermined for most libraries
  • The underlying linear algebra is static
aztecoo linear solver library
AztecOO Linear Solver Library
  • A C++ wrapper around Aztec library written in C.
  • Algorithms: GMRES, CG, CGS, BiCGSTAB, TFQMR.
  • Offers status testing capabilities.
  • Output verbosity level can be determined by user.
  • Uses Epetra to perform underlying vector space operations.
  • Interface to matrix-vector product is defined by the user through the Epetra_Operator.
arnoldi package arpack
ARnoldi PACKage (ARPACK)
  • Written in Fortran 77.
  • Algorithms: Implicitly Restarted Arnoldi/Lanczos
  • Static convergence tests.
  • Output formatting, verbosity level is determined by user.
  • Uses LAPACK/BLAS to perform underlying vector space operations.
  • Offers abstract interface to matrix-vector products through reverse communication.
slide6

Scalable Library for Eigenvalue Problem Computations( SLEPc )

  • Written in C (Hernández, Román, and Vidal, 2003).
  • Provides some basic eigensolvers as well as wrappers around:
    • ARPACK(Lehoucq, Maschhoff, Sorensen, and Yang, 1998)
    • BLZPACK (Marques, 1995)
    • PLANSO(Wu and Simon 1997)
    • TRLAN(Wu and Simon, 2001)
  • Native Algorithms: Power/Subspace Iteration, RQI, Arnoldi
  • Wrapped Algorithms: IRAM/IRLM (ARPACK), Block Lanczos (BLZPACK), and Lanczos (PLANSO / TRLAN)
  • Static convergence tests.
  • Uses PETSc to perform underlying vector space operations, matrix-vector products, and linear solves.
  • Allows the creation / registration of new matrix-vector products.
linear eigensolver software design
Linear/Eigensolver Software Design

Why not create a solver library that:

  • Provides an abstract interface to an operator-vector products, scaling, and preconditioning.
  • Allows the user to enlist any linear algebra package for the elementary vector space operations essential to the algorithm. (Epetra, PETSc, etc.)
  • Allows the user to define convergence of any algorithm (a.k.a. status testing).
  • Allows the user to determine the verbosity level, formatting, and processor for the output.
generic iterative method

Status Test

Linear/Eigen Problem

Generic Operator

Interface

Generic Linear

Algebra Interface

Output Manager

Generic Iterative Method

AlgorithmA( … )

while ( myStatusTest.CheckStatus(…)

== Unconverged )

% Compute operator-vector product

myProblem.ApplyOp(NOTRANS,v,w);

 = w.dot(v);

end

myOM.print(something);

return (solution);

end

benefits of generic programming
Benefits of Generic Programming
  • Generic programming techniques ease the implementation of complex algorithms.
  • Developing algorithms with generic programming techniques is easier on the developer, while still allowing them to build a flexible and powerful software package.
  • Generic programming techniques also allow the user to leverage their existing software investment.

Caveat: It’s not as easy as taking a piece of code and adding:

template <typename OT, typename ST>

More work has to be done to handle “numeric traits”

teuchos numeric traits
Teuchos Numeric Traits
  • OrdinalTraits
    • zero & one
    • int & long int
  • ScalarTraits
    • zero, one, magnitude type, absolute value, conjugate,

square root, random number generator, …

    • std::numeric_limits
    • float, double, complex<float>, and complex<double>
  • Arbitrary precision arithmetic
  • Templated BLAS/LAPACK wrappers, serial dense matrix/vector

Generic programming technique that uses templated interfaces to define the standard behavior of datatypes.

teuchos templated blas example
Teuchos Templated BLAS Example

template<typename OrdinalType, typename ScalarType>

ScalarType BLAS<OrdinalType, ScalarType>::NRM2( const OrdinalType n,

const ScalarType* x,

const OrdinalType incx) const

{

OrdinalType izero = OrdinalTraits<OrdinalType>::zero();

OrdinalType ione = OrdinalTraits<OrdinalType>::one();

ScalarType result = ScalarTraits<ScalarType>::zero();

OrdinalType i, ix = izero;

if ( n > izero )

{

// Set the initial index.

if (incx < izero) { ix = (-n+ione)*incx; }

for(i = izero; i < n; i++)

{

result += ScalarTraits<ScalarType>::conjugate(x[ix]) * x[ix];

ix += incx;

}

result = ScalarTraits<ScalarType>::squareroot(result);

}

return result;

} /* end NRM2 */

double DNRM2( const int n, const double* x, const int incx) const

{

int i, ix = 0;

double result 0.0;

if ( n > 0 )

{

// Set the initial index.

if (incx < 0) { ix = (-n+1)*incx; }

for(i = 0; i < n; ++i)

{

result += x[ix] * x[ix];

ix += incx;

}

result = std::sqrt(result);

}

return result;

} /* end NRM2 */

new packages belos and anasazi
New Packages: 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 (mini interface). Interfaces are derived from these base classes to operator-vector products, status tests, and any arbitrary linear algebra library.
  • Includes block linear solvers and eigensolvers.
  • Release Dates:
    • Internal (SNL): May 2004
    • External (Public): March 2005
why are block solvers useful
Why are Block Solvers Useful?
  • Block Solvers ( in general ):
    • Achieve better performance for operator-vector products.
  • Block Eigensolvers ( Op(A)X = X ):
  • Block Linear Solvers ( Op(A)X = B ):
  • Reliably determine multiple and/or clustered eigenvalues.
  • Example applications: Modal analysis, stability analysis,

bifurcation analysis (LOCA)

  • Useful for when multiple solutions are required for the same system of equations.
  • Example applications:
    • Perturbation analysis
    • Optimization problems
    • Single right-hand sides where A has a handful of small eigenvalues
    • Inner-iteration of block eigensolvers
belos anasazi status trilinos release 4 0
Belos/Anasazi Status(Trilinos Release 4.0)
  • Belos:
    • Solvers: BlockCG and BlockGMRES
    • Linear algebra adapters for Epetra and NOX/LOCA
    • Belos::EpetraOperator allows for integration into other codes
    • Epetra interface accepts Epetra_Operators, so can be used with ML, AztecOO, Ifpack, Belos, etc…
    • Block size is independent of number of right-hand sides
  • Anasazi:
    • Solvers: Block Krylov-Schur Method and Modal Analysis Solvers (Davidson, Jacobi-Davidson, LOBPCG, etc…)
    • Linear algebra adapters for Epetra and NOX/LOCA
    • Can solve standard and generalized eigenproblems
    • Epetra interface accepts Epetra_Operators, so can be used with Belos, AztecOO, etc…
    • Block size is independent of number of requested eigenvalues
belos anasazi design
Belos/Anasazi Design
  • LinearSolver / Eigensolver Class
  • LinearProblem / Eigenproblem Class
  • Sorting and Orthogonalization Class
  • StatusTest Class
  • OutputManager Class
  • Configurable via Teuchos::ParameterList
  • Interface current Belos/Anasazi abstract linear algebra classes with TSFCore abstract linear algebra classes.
belos anasazi linear algebra interface mini interface
Belos/Anasazi Linear Algebra Interface ( Mini-interface )
  • (Belos/Anasazi)::MultiVec
    • Abstract interface to define the linear algebra required by most iterative linear solvers / eigensolvers :
      • clone methods
      • dot product
      • update methods
      • norms
      • initialize / randomize
  • (Belos/Anasazi)::Operator
    • Abstract interface to define the linear operator interface required by most iterative linear solvers / eigensolvers :
      • apply methods
belos anasazi iterativesolver interface
Belos/Anasazi IterativeSolver Interface
  • Provides an abstract interface to Belos/Anasazi basic iterations.
  • (Belos/Anasazi)::IterativeSolver
    • number of iterations / restarts
    • native residuals
    • linear problem / eigenproblem
    • initialize/finalize
    • iterate
  • Specialization for solvers that can provide more information
    • Block Krylov-Schur -> Ritz values / residuals
belos anasazi linear eigenproblem interface
Belos/Anasazi Linear/Eigenproblem Interface
  • Provides an interface between the basic iterations and linear/eigenproblem to be solved.
  • Belos::LinearProblem
    • Allows preconditioning, scaling, etc. to be removed from the algorithm, performs composition apply.
    • Preconditioner/scaling, operator, LHS/RHS, StatusTest
    • Performs bookkeeping for solver
  • Anasazi::Eigenproblem
    • Allows spectral transformations to be removed from the algorithm.
    • Differentiates between standard and generalized eigenproblems
    • Intial vector, stiffness/mass matrix, operator, eigenvalues, eigenvectors
belos anasazi statustest interface
Belos/Anasazi StatusTest Interface
  • Similar to NOX / AztecOO.
  • Allows user to decide when solver is finished.
  • Multiple criterion can be logically connected.
  • Abstract base class methods:
    • StatusType CheckStatus( IterativeSolver<Scalar>* solver );
    • StatusType GetStatus() const;
    • void Reset();
    • ostream& Print( ostream& os ) const;
  • Derived classes
    • Maximum iterations/restarts,
    • Residual norm,
    • Combinations {AND,OR}, …
belos anasazi output manager interface
Belos/Anasazi Output Manager Interface
  • Allows user to control which processor will handle output from the solver and the verbosity.
  • Default is lowest verbosity, outputting on proc 0.
  • Methods:
    • Get/Set Methods:
      • void SetOStream( ostream& os );
      • void SetVerbosity( int vbLevel );
      • ostream& GetOStream();
    • Query Methods:
      • bool doOutput( int vbLevel );
anasazi sort manager interface
Anasazi Sort Manager Interface
  • Abstract interface for managing the sorting of the eigenvalues computed by the eigensolver
  • Important tool to complement spectral transformations
  • Only two methods:
    • ReturnType sort(IterativeSolver<Scalar>* solver, int n, TYPE *evals, int *perm=0);
    • ReturnType sort(IterativeSolver<Scalar>* solver, int n, TYPE *r_evals, TYPE *i_evals, int *perm=0);
  • Basic sorting types are provided
    • largest/smallest magnitude, largest/smallest real part, largest/smallest imaginary part
how do the interfaces interoperate anasazi
How do the Interfaces Interoperate?(Anasazi)

Anasazi::IterativeSolver

Anasazi::SortManager

Anasazi::OutputManager

Anasazi::Eigenproblem

Anasazi::MultiVec

Anasazi::Operator

Anasazi::StatusTest

projected belos anasazi status trilinos release 5 0
Projected Belos/Anasazi Status(Trilinos Release 5.0)
  • Belos:
    • Solvers: CG, BlockCG, BlockGMRES, and BlockFGMRES
    • Interfaces are defined for Epetra, NOX/LOCA, and TSFCore
    • Configurable via a Teuchos::ParameterList
  • Anasazi:
    • Solvers: Block Krylov-Schur Method and Modal Analysis Solvers (Davidson, Jacobi-Davidson, LOBPCG, etc…)
    • Interfaces are defined for Epetra, NOX/LOCA, and TSFCore
    • Abstract interface to sorting and orthogonalization (???)
    • Configurable via a Teuchos::ParameterList
future anasazi belos development
Future Anasazi/Belos Development
  • Belos
    • Algorithms: FCG and BiCGSTAB
    • GMRES not necessary
    • Abstract orthogonalization class, user-determined Krylov subspace breakdown
    • Higher-level robust managers to orchestrate solvers
  • Anasazi
    • Algorithms: Restarted Preconditioned Eigensolver (RPE), Generalized Davidson
    • Abstract orthogonalization class, user-determined Krylov subspace breakdown
    • Higher-level robust managers to orchestrate solvers
concluding remarks audience participation
Concluding Remarks(audience participation)
  • Are there any capabilities you may need from a linear/eigensolver that have not been addressed?
  • Do you have any questions, comments, or suggestions about the proposed user interface?
  • Are there any other solvers that you think would be useful?