NOX: An Object-Oriented Nonlinear Solver Package Roger Pawlowski , Tamara Kolda, Russell Hooper,

1. NOX: An Object-Oriented Nonlinear Solver Package Roger Pawlowski, Tamara Kolda, Russell Hooper, and John Shadid Sandia National Laboratories Albuquerque, NM Trilinos User Group Meeting October 16th, 2003 Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company,for the United States Department of Energy’s National Nuclear Security Administration under contract DE-AC04-94AL85000.

2. Background • NOX is a Trilinos solver package built to: • enable robust and efficient solutions to systems of nonlinear equations. • provide a broad array of nonlinear algorithms through a single interface. • rapid deployment of new solver technology into ASCI codes. • eliminate redundancies in ASCI code development. • Started ~2.5 years ago. • Funding: ASCI Algorithms and LDRD (2.5 FTE). • Consists of two libraries: • NOX – nonlinear algorithms • LOCA v2.0 – continuation and bifurcation algorithms (Eric Phipps) • Website: software.sandia.gov/nox

3. In-NOX-ulated Applications Primary NOX Team: Tammy Kolda, Roger Pawlowski, Russ Hooper

4. Nonlinear Equations • Given • find for which • F is a system of n nonlinear equations with n unknowns. • x is the unknown or solution vector of size n.

5. Tensor Method Newton’s Method Broyden’s Method Line Search Interval HalvingQuadratic Cubic More’-Thuente Homotopy Artificial Parameter Continuation Natural Parameter Continuation Trust Region Dogleg Inexact Dogleg Nonlinear Solution Algorithms Globalizations Iterative Linear Solvers, Adaptive Forcing Terms User Norms and Merit Functions

6. Building Blocks of NOX Calculating the Direction Iterate Control (Solver) Damping or Line Search Stopping Criteria (Status Test) Example: Newton’s Method for F (x)= 0 • Choose an initial guess x0 • For k = 0,1,2,... • Compute Fk = F (xk) • Compute Jk where (Jk )ij = ¶F i(xk)/¶x j • Let dk = -Jk-1Fk • (Optional) Let lk be a calculated step length • Set xk+1= xk + lkdk • Test for Convergence or Failure

7. NOX Solvers(Iterate Control) Derived from NOX::Solver::Generic bool  reset(NOX::Abstract::Group& grp, NOX::StatusTest::Generic& tests, NOX::Parameter::List& params) bool  reset(NOX::Abstract::Group& grp, NOX::StatusTest::Generic& tests) NOX::StatusTest::StatusType getStatus() NOX::StatusTest::StatusType iterate() NOX::StatusTest::StatusType solve() const NOX::Abstract::Group&  getSolutionGroup() const const NOX::Abstract::Group&  getPreviousSolutionGroup() const int getNumIterations() const const NOX::Parameter::List&  getParameterList() const • LineSearchBased • Compute Direction • Compute Step Length • Scale direction by step length • Update solution • Check Convergence • TrustRegionBased (dogleg) • Computes Newton Direction • Computes Cauchy Direction • Adjust trust region radius • Update solution • Check Convergence • InexactTrustRegionBased (Prerelease) • TensorBased(Prerelease) • UserDefined Solver is determined by passing a parameter list to a handle for the class: NOX::Solver::Manager solver(group, StatusTests, ParameterList); int status = solver.solve();

8. NOX Directions • Steepest Descent • Newton • Requires a linear solve • Adaptive forcing terms • Broyden • Modified Newton with rank-1 updates to the most recently computed Jacobian • Nonlinear-CG (Prerelease) • Tensor (Prerelease) Direction is determined and constructed in a handle called a “Manager”: NOX::Direction::Manager dir(PrintParams, DirectionParams); User Defined Directions Derived from NOX::Direction::Generic Passed into the solver via parameter list using a direction template ctor. User Defined Merit Function Derived from NOX::Parameter::MeritFunction Passed into the solver via parameter list using the NOX::Parameter::Arbitrary entry

9. NOX Line Searches Line Searches are determined and constructed in a handle called a “Manager”: NOX::LineSearch::Manager ls(PrintParams, LineSearchParams); • Full Step • Constant step size. • Defaults to 1.0. • Backtrack • interval halving • Polynomial Interpolation • Quadratic and Cubic interpolation • Requires a “Sufficient Decrease” Condition • More’-Thuente • Polynomial Interpolation • Sufficient Decrease and Curvature Conditions • Nonlinear-CG – (Prerelease) • Curvilinear (Prerelease) User Defined Line Searches Derived from NOX::LineSearch::Generic Passed into the solver via parameter list using a template constructor. User Defined Merit Function Derived from NOX::Parameter::MeritFunction Passed into the solver via parameter list using the NOX::Parameter::Arbitrary entry User Defined Norms Derived from NOX::Parameter::UserNorm Passed into the solver via parameter list using the NOX::Parameter::Arbitrary entry

10. Newton Quadratic Line Search Tensor Method Trust-Region Method Newly-Deployed Broyden Method Mix-n-Match Solver Algorithms • Higly versatile code environment • Solver • Line Search Based • Trust Region Based • Tensor Based • Direction • Newton • Broyden • Steepest Descent • Tensor • User-Defined • Line Search / Damping • Full Step • Backtrack • Polynomial/Quadratic • More’-Thuente • Curvilinear • User-Defined

11. Stopping Criteria (StatusTests) • Highly Flexible Design: Users build a convergence test hierarchy and registers it with the solver (via solver constructor or reset method). • Norm F: {Inf, One, Two} {absolute, relative} • Norm Update DX: {Inf, One, Two} • Norm Weighted Root Mean Square (WRMS) • Max Iterations – Failure test if solvers reaches max iters • FiniteValue – Failure test that checks for NaN and Inf on • Stagnation – Failure test that triggers if the convergence rate fails a tolerance check for n consecutive iterations. • Combination: {AND, OR} • Users Designed: Derive from NOX::StatusTest::Generic Return Types Unconverged Converged* Failed* Unevaluated

12. FiniteValue: finiteValueTest Combo(OR) allTests MaxIters: maxItersTest normFTest Combo(AND): convergedTest normWRMSTest NOX::StatusTest::NormF normFTest(); NOX::StatusTest::NormWRMS normWRMSTest(); NOX::StatusTest::Combo convergedTest(NOX::StatusTest::Combo::AND); convergedTest.addStatusTest(normFTest); convergedTest.addStatusTest(normWRMSTest); NOX::StatusTest::FiniteValue finiteValueTest; NOX::StatusTest::MaxIters maxItersTest(200); NOX::StatusTest::Combo allTests(NOX::StatusTest::Combo::OR); allTests.addStatusTest(finiteValueTest); allTests.addStatusTest(maxItersTest); allTests.addStatusTest(convergedTest); Building a Status Test • Converge if both: • Fail if value of becomes Nan or Inf • Fail if we reach maximum iterations

13. Status Tests Continued -- Status Test Results --**...........OR Combination -> **...........AND Combination -> **...........F-Norm = 5.907e-01 < 1.000e-08 (Length-Scaled Two-Norm, Absolute Tolerance) **...........WRMS-Norm = 4.794e+01 < 1 (Min Step Size: 1.000e+00 >= 1) (Max Lin Solv Tol: 1.314e-15 < 0.5) **...........Finite Number Check (Two-Norm F) = Finite **...........Number of Iterations = 2 < 200 -- Final Status Test Results --Converged....OR Combination -> Converged....AND Combination -> Converged....F-Norm = 3.567e-13 < 1.000e-08 (Length-Scaled Two-Norm, Absolute Tolerance) Converged....WRMS-Norm = 1.724e-03 < 1 (Min Step Size: 1.000e+00 >= 1) (Max Lin Solv Tol: 4.951e-14 < 0.5) ??...........Finite Number Check (Two-Norm F) = Unknown ??...........Number of Iterations = -1 < 200 User Defined are Derived from NOX::StatusTest::Generic NOX::StatusTest::StatusType checkStatus(const NOX::Solver::Generic &problem) NOX::StatusTest::StatusType  checkStatusEfficiently(const NOX::Solver::Generic &problem, NOX::StatusTest::CheckType checkType) NOX::StatusTest::StatusType getStatus() const ostream& print(ostream &stream, int indent=0) const

14. Solvers- Line Search - Trust Region Directions- e.g., Newton Line Searches - e.g., Polynomial Status Tests - e.g., Norm F AbstractLayer Abstract Vector & Abstract Group NOX Framework SolverLayer • Don’t need to directly access the vector or matrix entries, only manipulate the objects. • NOX uses an abstract interface to manipulate linear algebra objects. • Isolate the Solver layer from the linear algebra implementations used by the application. • This approach means that NOX does NOT rely on any specific linear algebra format. • Allows the apps to tailor the linear algebra to their own needs! • Serial or Parallel • Any Storage format: LAPACK, PETSc, Epetra, User Defined

15. Initialization x = y x = |y| xi = 1/yifor i = 1 to n xi = g for i = 1 to n Length of Vector Scaling x = ax xi = xiyifor i = 1 to n Update x = aa + gx, x = aa + bb + gx Norm kxk1, kxk2 , kxk1 kxkw (weighted norm) Dot x¢y Clone (create a copy) y = x NOX::Abstract::Vector The solver is not allowed nor does it need explicit access to the vector – just the ability to manipulate it.

16. x = Iterate Initialize Update Access F = F (x) Compute Access J = Jacobian of F at x Compute Apply / Apply Transpose Apply Inverse Apply Preconditioning n = Newton Vector Compute n = -J-1F to specified tolerance Access g = Gradient of kF (xk)k2 g = JTF Access Clone NOX::Abstract::Group The linear solver and application interface are combined into the Group. Vectors x, F, n, g are accessed as NOX::Abstract::Vector’s. Matrix J is never directly accessed!

17. Code Demonstration

18. Solvers- Line Search - Trust Region Directions- e.g., Newton Line Searches - e.g., Polynomial Status Tests - e.g., Norm F AbstractLayer Abstract Vector & Abstract Group Implementations- EPetra- PETSc - LAPACK- USER DEFINED Linear Algebra Interface EPetra Dependent Features- Matrix-Free Newton-Krylov- Preconditioning- Graph Coloring / Finite Diff. User Interface- Compute F- Compute Jacobian- Compute Preconditioner ApplicationInterface Layer NOX Framework SolverLayer

19. Linear Algebra Support Features • NOX’s support libraries define a concrete implementation of the Abstract Vector and Group. • We only require the user to implement a minimal interface derived from NOX::<Package>::Interface. bool computeF(const Epetra_Vector &x, Epetra_Vector &f, FillType flag=F) bool computeJacobian(const Epetra_Vector &x, Epetra_Operator &Jac) bool computePrecMatrix(const Epetra_Vector &x, Epetra_RowMatrix &M) bool computePreconditioner(const Epetra_Vector &x, Epetra_Operator &M) Group(NOX::Parameter::List& printingParams, NOX::Parameter::List& linearSolverParams, NOX::Epetra::Interface& i, NOX::Epetra::Vector& x, Epetra_Operator& J) Group(NOX::Parameter::List& printingParams, NOX::Parameter::List& linearSolverParams, NOX::Epetra::Interface& i, NOX::Epetra::Vector& x, Epetra_Operator& J, Epetra_Operator& M)

20. The Epetra “Goodies” • Matrix-Free Newton-Krylov Operator • Derived from Epetra_Operator • Can be used to estimate Jacobian action on a vector • NOX::Epetra::MatrixFree • Finite Difference Jacobian • Derived from an Epetra_RowMatrix • Can be used as a preconditioner matrix • NOX::Epetra::FiniteDifference • Graph Colored Finite Difference Jacobian • Derived from NOX::Epetra::FiniteDifference • Fast Jacobian fills – need connectivity/coloring graph • (NOX::Epetra::FiniteDifferenceColoring) • Full interface to AztecOO using NOX parameter list • Preconditioners: internal AztecOO, Ifpack, User defined • Scaling object

21. Parameter Continuation Solve a series of problems that progressively become the system of interest based on a parameter, a: Natural: FH (x) = F(x,a) Artificial: FH (x) = aF(x) + (1 – a)G(x) F(x) G(x) easy 0 1 a Homotopy Algorithms • LOCA: Library of Continuation Algorithms • Bifurcation and stability analysis package. • Tightly coupled to NOX to reuse application interface. • Provided as a single package. • Stepper object provides a robust and intelligent step control for parameter continuation. • Eric Phipps

22. NOX Specific Configure Options • Compiling NOX library in Trilinos: • --enable-nox • Compiling prerelease code: • --enable-prerelease • Using nox built-in linear algebra support: • LAPACK: --enable-nox-lapack • --enable-nox-lapack-examples • Epetra: --enable-nox-epetra • --enable-nox-epetra-examples • PETSc: --enable-nox-petsc • --enable-nox-petsc-examples • Compiling the test suite: • --enable-tests • --enable-nox-tests • Compiling LOCA library in Trilinos/NOX: • --enable-loca

23. Future Work • Test Suite Development – Darin Diachin • Re-work of the NOX::Epetra support library for efficiency - Pawlowski • TSF support library – Pawlowski, Hooper • Tensor algorithm development – Brett Bader, Tammy Kolda • Broyden algorithm development – Kolda, Pawlowski • Multi-Physics Solvers – Strong code coupling

24. Multi-Physics Solvers • MultiPhysics – Bill Spotz and Alfred Lorber • rapid prototyping environment to explore solver coupling algorithms. • PyNOX - Python interface to call NOX solver. • SIERRA – Russell Hooper, John Shadid and Roger Pawlowski • Implementation of a Jacobian-Free coupling algorithm in SIERRA. Brusselator Example in NOX (Russell Hooper) Loose Coupling Jacobian-Free Coupling Full Coupling with preconditioner: Analytic OR FD Coloring (Epetra)

25. Summary • Key Features • Object-oriented C++ library • Large number of cutting-edge algorithms • Can use any linear algebra implementation • Easy interface • User Flexibility: Add solvers, directions, line searches, convergence tests • NOX Support libraries • Fast integration – reduced interface • parallel iterative solvers (AztecOO, PETSc) • preconditioners (Ifpack, ML, Aztec, PETSc) • Matrix-Free Newton-Krylov • Matrix Estimation: Finite Difference/Coloring

26. Sandia Tammy Kolda * Roger Pawlowski * Russ Hooper * John Shadid Todd Coffey * Andy Salinger (LOCA) * Eric Phipps (LOCA) * Bill Spotz * Mike Heroux (Trilinos) Scott Hutchinson (Xyce) Rob Hoekstra (Xyce) Eric Keiter (Xyce) Alan Williams (SIERRA) Tom Smith (Premo) Alfred Lorber (Premo) Tom Brunner (Alegra) Academics Homer Walker, WPI Joe Simonis, WPI Tim Kelley, NCSU Bobby Schnabel, UCB Richard Byrd, UCB Ryan McKenzie, Kentucky Craig Douglas, Kentucky NOX Contributors (Ideas, Code, and Testing) * NOX Developer NOX Information Tammy Kolda: tgkolda@sandia.gov Roger Pawlowski: rppawlo@sandia.gov Download: http://software.sandia.gov/nox