Créer une présentation
Télécharger la présentation

Download

Download Presentation

An Overview of Belos and Anasazi October 16 th , 11-12am Heidi Thornquist Teri Barth Rich Lehoucq Mike Heroux Computatio

302 Vues
Download Presentation

Télécharger la présentation
## An Overview of Belos and Anasazi October 16 th , 11-12am Heidi Thornquist Teri Barth Rich Lehoucq Mike Heroux Computatio

- - - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - - -

**An Overview of Belosand AnasaziOctober 16th, 11-12amHeidi**ThornquistTeri BarthRich LehoucqMike HerouxComputational Mathematics and Algorithms 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**• Background • AztecOO and ARPACK • Solver design requirements • Belos and Anasazi • Current Status • Current Development • Teuchos • Arbitrary precision BLAS • Concluding remarks**Background**• Several iterative linear solver / eigensolver libraries exist: • PETSc, SLAP, LINSOL, Aztec(OO), … • LOBPCG (hypre), ARPACK, … • None are written in C++ • None of the linear solver libraries can efficiently deal with multiple right-hand sides • Current eigensolver libraries offer one algorithm for solving symmetric/non-symmetric eigenproblems**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 EpetraOperator.**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.**Generic programming**Some design requirements for Linear Solvers /Eigensolvers • Abstract interface to an operator-vector product (Op(A)v) and preconditioning (if necessary). • Ability to allow the user to enlist any linear algebra package for the elementary vector space operations essential to the algorithm • Ability to allow the user to define convergence of any algorithm (a.k.a. status testing). • Ability to allow the user to determine the verbosity level, formatting, and processor for the output.**Status Test**Generic Problem Interface Generic Linear Algebra Interface IO Manager Generic Algorithm AlgorithmA( … ) while ( myStatusTest.GetStatus(…) != converged ) … … % Compute operator-vector product myProblem.ApplyOp(NOTRANS,v,w); = w.dot(v); … end myIO.print(something); return (solution); end**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.**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. Interfaces are derived from these base classes to operator-vector products, status tests, and any arbitrary linear algebra library. • Includes block linear solvers (GMRES, CG, etc.) and eigensolvers (IRAM, etc.). • Release: Spring/Summer 2004**Why are Block Solvers Useful?(from an applications point of**view) • Block Eigensolvers ( Op(A)X = X ): • Block Linear Solvers ( Op(A)X = B ): • Reliably determine multiple and/or clustered eigenvalues. • Example applications: • 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**Why are Block Solvers Useful?(from a computational point of**view) • Remember: These are iterative solvers! • Solutions for block solvers are obtained through operator-multivector computations (Op(A)V). • However, single-vector iterative solvers use operator-vector computations (Op(A)v). • What’s the difference? BLAS Performance • Op(A)v is Level 1 • Op(A)V is Level 2**Current Status(Trilinos Release 3.2d)**• Belos: • Solvers: BlockCG and BlockGMRES • Interfaces are defined for Epetra and NOX/LOCA • BelosEpetraOperator allows for integration into other codes • Epetra interface accepts EpetraOperators, so can be used with ML, AztecOO, Ifpack, Belos, etc… • Block size is independent of number of right-hand sides • Anasazi: • Solver: Implicitly Restarted Arnoldi • Interfaces are defined for Epetra and NOX/LOCA • Can solve standard and generalized eigenproblems • Epetra interface accepts EpetraOperators, so can be used with Belos, AztecOO, etc… • Block size is independent of number of requested eigenvalues**Current Belos/Anasazi Development**• LinearSolver / Eigensolver Class • LinearProblem / Eigenproblem Class • StatusTest Class • IOManager Class • Incorporate parameter lists ( developed in Teuchos ). • Replace current Belos/Anasazi abstract linear algebra classes with TSFCore abstract linear algebra classes.**Current Belos/Anasazi DevelopmentLinearSolver / Eigensolver**Class • Provide an abstract interface to Belos/Anasazi solvers. • Methods: • Set Methods: SetProblem, SetParameters, SetStatusTester • Get Methods: GetProblem, GetParameters, GetStatusTester, currIteration • Solve Methods: Solve, Iterate • Auxilliary Method(s): Reset**Belos**Allows preconditioning, scaling, etc. to be removed from the algorithm. Methods: Set Methods: SetOperator, SetLHS, SetRHS, SetL/RScale, SetL/RPrecond Get Methods: GetOperator, GetLHS, GetRHS, GetL/RPrecond Query Methods: IsL/RPrecond Apply Methods: ApplyOperator Anasazi Differentiates between standard and general eigenvalue problems. Allows spectral transformations to be removed from the algorithm. Methods: Set Methods: SetInitVec, SetOperator, SetMatrixA, SetMatrixB Get Methods: GetInitVec, GetOperator, GetMatrixA, GetMatrixB, GetEvals, GetEvecs, GetResids Apply Methods: ApplyOperator, ApplyMatrixA, ApplyMatrixB IP/Norm Methods: AInProd, BInProd, BMvNorm Current Belos/Anasazi DevelopmentLinearProblem / Eigenproblem Class**Current Belos/Anasazi DevelopmentStatusTest Class**• Similar to NOX / AztecOO. • Allows user to decide when solver is finished. • Multiple criterion can be logically connected. • Abstract base class methods: • StatusType CheckStatus( int currIteration, LinearProblem<Scalar>& myProblem ); • StatusType CheckStatus( int currIteration, Eigenproblem<Scalar>& myProblem ); • StatusType GetStatus() const; • ostream& Print( ostream& os ) const; • Derived classes • Maximum iterations/restarts, residual norm, combinations {AND,OR}, …**Current Belos/Anasazi DevelopmentIOManager Class**• Allows user to control I/O processor, precision, and verbosity. • Methods: • Set Methods: • void reset(Teuchos::ParameterList& myIOParams ); • Query Methods: • bool isPrintProc() const; • bool isPrintLevel() const; • bool isPrintProcLevel() const; • Output Methods: overload “<<“ operator using given precision.**Belos**CG/FCG GMRES/FGMRES BiCGSTAB Anasazi LOBPCG Restarted Preconditioned Eigensolver (RPE) Generalized Davidson Future Belos/Anasazi Algorithms**Belos**Preprocessing capabilities (CG) User determined Krylov-subspace breakdown tolerances Anasazi Preprocessing capabilties User determined sorting User determined Krylov-subspace breakdown tolerances Future Belos/Anasazi Features**Belos / Anasazi Development(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?**Teuchos**• Optional, lightweight utility package • Takes advantage of advanced features of C++: • Templates • Standard Template Library (STL) • Planned features: • Parameter list of arbitrary types (NOX) • Memory management package and error-handler (TSFCoreUtils) • Templated interfaces to BLAS, LAPACK, dense matrix (TPetra) • Command line and XML parsers, random number generators, and timers (TSF)**BLAS**• Basic Linear Algebra Subroutines • Robust, high-performance dense vector-vector, matrix-vector, and matrix-matrix routines • Written in Fortran • Fast and efficient • Fortran is highly portable; however, the interface between it and C++ is platform-dependent • Templated C++ wrappers for the BLAS Fortran routines are created through autoconf. • Two main benefits of using a templated interface • Template Specialization: template specialization can be used to support the four native datatypes of existing Fortran BLAS code • General Implementation: a generic C++ implementation can be created to handle many arbitrary datatypes**Arbitrary Datatypes**• Any datatype that defines zero, one, addition, subtraction, and multiplication can use nearly all BLAS functions • Some require square root or division • Example: ARPREC • arbitrary precision library that uses arrays of 64-bit floating-point numbers to represent high-precision floating-point numbers. • ARPREC values behave just like any other floating-point datatype, except the maximum working precision (in decimal digits) must be specified before any calculations are done mp::mp_init(200); • Illustrate the use of ARPREC with an example using Hilbert matrices.**Hilbert Matrices**• A Hilbert matrix HN is a square N-by-N matrix such that: • Notoriously ill-conditioned • κ(H3) 524 • κ(H5) 476610 • κ(H10) 1.6025 x 1013 • κ(H20) 7.8413 x 1017 • κ(H100) 1.7232 x 1020 • Theoretically, Hilbert matrices are SPD**Cholesky Factorization and Hilbert Matrices**• In practice, floating-point rounding error causes higher-order Hilbert matrices to no longer be positive definite • Cholesky factorization is often done to determine whether or not a given matrix is positive definite • As such, a Cholesky factorization will often fail on higher-order Hilbert matrices, even though such matrices “should” be positive definite**Cholesky Factorization and Hilbert Matrices**• In addition to providing a higher limit on N, using arbitrary precision reduces error. • By solving the linear system: where the vector b is constructed such that the true solution x is equal to a 1-vector we can look at the forward error.**Templated BLAS Conclusions**• Adding a templated interface to the BLAS allows the use of arbitrary datatypes • This is useful in particular for the inclusion of arbitrary-precision arithmetic • Such arithmetic may be expensive, but the increased precision may be required for ill-conditioned problems • Arbitrary-precision arithmetic can be used in only the places where it is needed most in order to retain speed • Templated tools libraries are essential for the development of templated solver libraries.