llnl / snls

Small Non-Linear Solver

	SNLS -- small non-linear solver library

 ________  ________   ___       ________      
|\   ____\|\   ___  \|\  \     |\   ____\     
\ \  \___|\ \  \\ \  \ \  \    \ \  \___|_    
 \ \_____  \ \  \\ \  \ \  \    \ \_____  \   
  \|____|\  \ \  \\ \  \ \  \____\|____|\  \  
    ____\_\  \ \__\\ \__\ \_______\____\_\  \ 
   |\_________\|__| \|__|\|_______|\_________\
   \|_________|                   \|_________|

BACKGROUND

SNLS (small non-linear solver) is a C++ library for solving non-linear systems of equations. The emphasis is on small systems of equations for which direct factorization of the Jacobian is appropriate. The systems of equations may be numerically ill-conditioned. Implementation is amenable to GPU usage, with the idea that class instantiation would occur within device functions for non-batch type solvers.

Examples of work that has made use of SNLS or substantially equivalent algorithms:

• Journal publication on use with a porosity mechanics type model
• Journal publication on use with a crystal mechanics type model

Solvers currently in the library:

• SNLSTrDlDenseG Dogleg approximation to the trust-region sub-problem for multi-dimensional nonlinear systems of equations. This method reduces to a Newton-Raphson approach near the solution. These methods are sometimes described as "restricted step" approaches instead of trust-region approaches. See, for example, Practical Methods of Optimization. As compared to general trust-region approaches in scalar minimization, the approach here for solving non-linear systems amounts to assuming that the Hessian matrix can be approximated using information from the Jacobian of the system.
• SNLSHybrdTrDLDenseG A hybrid trust region type solver, dogleg approximation for dense general Jacobian matrix that makes use of a rank-1 update of the jacobian using QR factorization. The method is inspired by SNLS current trust region dogleg solver and Powell's original hybrid method for nonlinear equations, and MINPACK's modified version of it. Powell's original hybrid method can be found at: M. J. D. Powell, "A hybrid method for nonlinear equations", in Numerical methods for nonlinear algebraic equations. MINPACK's user guide is found at https://doi.org/10.2172/6997568 . As compared to general trust-region approaches in scalar minimization, the approach here for solving non-linear systems amounts to assuming that the Hessian matrix can be approximated using information from the Jacobian of the system.
• SNLSTrDlDenseG_Batch a batch version of SNLSTrDlDenseG that is only available if the library is compiled with the -DUSE_BATCH_SOLVERS=On cmake option. It makes use of the RAJA Performance Suite to abstract away the complexity of memory management and execution strategies of forall loops. The test/SNLS_batch_testdriver.cc provides an example of how this solver, the memory manager, and forall abstraction layer can be used to solve a modified Broyden Tridiagonal Problem.
• NewtonBB Simple 1D Newton solver with a fallback to bisection. If the zero of the function can be bounded then this solver can return a result even if the function is badly behaved.

BUILDING

The build system is cmake-based.

Dependencies:

• blt -- required
  • https://github.com/LLNL/blt
  • in cmake invocation specify location with -DBLT_SOURCE_DIR

TESTING

Run cmake with -DENABLE_TESTS=ON and do make test

DEVELOPMENT

The develop branch is the main development branch for snls. Changes to develop are by pull request.

AUTHORS

The principal devleoper of SNLS is Nathan Barton, [email protected]. Brett Wayne and Robert Carson have also made contributions.

CITATION

SNLS may be cited using the following bibtex entry:

@Misc{snls,
title = {SNLS},
author = {Wayne, Brett M. and Barton, Nathan R. and USDOE National Nuclear Security Administration},
abstractNote = {{SNLS} (small non-linear solver) is a C++ library for solving non-linear systems. The emphasis is on small systems of equations for which direct factorization of the Jacobian is appropriate. The systems of equations may be numerically ill-conditioned. Implementation is amenable to {GPU} usage, with the idea that class instantiation would occur within device functions.},
doi = {10.11578/dc.20181217.9},
year = 2018,
month = {9},
url = {https://github.com/LLNL/SNLS},
annote =    {
   https://www.osti.gov//servlets/purl/1487196
   https://www.osti.gov/biblio/1487196
}
}

LICENSE

License is under the BSD-3-Clause license. See LICENSE file for details. And see also the NOTICE file.

SPDX-License-Identifier: BSD-3-Clause

LLNL-CODE-761139

snls's Issues

Update README to mention new solver

Realized I missed this as part of #6 ... As the HIP stuff gets closer to being ready (probably within this month), I can just include the change in there rather than opening up a new PR.

Missing getNJEvals function

SNLS/src/SNLS_TrDLDenseG.h

Line 82 in ca6f8e9

__snls_hdev__ int getNFEvals() const { return(_fevals ); };

The Hybrid solver has an additional member function for querying the number of Jacobian evaluations.

Add HIP backened to run on AMD hardware

A few changes are needed in order to make use of the HIP compiler suite. I'll work on this on a separate branch and then make the appropriate PR to bring this back in.

Potential further abstractions for working arrays / raja views

So on occasions, a regular RAJA view associated with a chai working array may not be the most convenient as we might like to also use that data on a different memory execution space for some purpose (debugging info, logging info, or something else like a calculation can only take on host for example). I'm contemplating a more generic class that can hold RAJA views for both the host and device class and then swap over to the appropriate one automatically based on the RAJA::plugins. I would imagine being able to do this in a manner similar to how CHAI's ManagedArrays work in automatically swapping things and not really too complicated to implement. Mainly I see this being useful in codes that make use of SNLS's device forall and memory manager abstraction layers. Also, I do see some areas it could be useful with the batch stuff in terms of debugging, or in the crazy edge case that someone decides to swap the execution strategy from GPU to CPU or GPU to CPU after the batch solver object has been created...

Potential Future Solvers/Algorithms

I'm just going to keep a list here of solvers/algorithms that might be interesting to implement in SNLS. I'll keep adding to this list as I find more that might be of interest.

I cam across this one today: https://arxiv.org/abs/2112.02089 which could be seen as a variation of a Levenberg-Marquardt type solvers. I was just looking through it (alg 2 mainly), and it seems like it would fairly simple to implement given it's not too far off from what we already do. If the convergence properties hold for our type of problems this could be a very interesting one to use especially for our very stiff equations.

Another one that might be of interest would be potentially making use of something like the double-dogleg in our solvers. You can find a description of it at section 6.4 of https://doi.org/10.1137/1.9781611971200.ch6 or https://doi.org/10.1007/BF00932218 . Although, it seems that this is largely the same as the single dogleg with the caveat that earlier on it would be a bit more biased towards the newton step rather than the cauchy point. However, I'm not necessarily sure it would be worth adding unless it drastically seems to help our solvers, which I don't believe it will.

Examine new solver exit conditions

The current solver just checks that the norm of the residual vector in the solve is less than some tolerance. I would say a more appropriate approach would be something similar to what's done in CG or other nonlinear solvers.

residual_norm_0 = ||b - Ax||
exit_tolerance = max(absolute_tolerance, relative_tolerance * residual_norm_0)

Additional methods related to this are outlined in https://doi.org/10.1137/0917003 which there's an implementation within MFEM's NewtonSolver https://github.com/mfem/mfem/blob/master/linalg/solvers.hpp#L460-L545

I noticed this becoming an issue in some tightly coupled PDE systems I'm solving for where the solver had clearly reached a steady state residual but the set tolerance was too high. I've had good experience with the 1st option, and it should be simple enough to implement within here. Although, we might want to think about the 2nd option as well. Since, it could allow for a more robust solver for certain sets of PDEs.

Add new tests for our N-dimension solvers

For our N-dimension solvers, we should look at maybe testing our solvers out with some of the test suites given here: https://github.com/devernay/cminpack/tree/master/examples and here: https://people.sc.fsu.edu/~jburkardt/f_src/test_nonlin/test_nonlin.html . These tests are from the MINPACK test suite for the nonlinear solvers and should at least provide confidence in any new N-dimension solvers that we might add. I only mention this as I've found that the Broyden test didn't catch an error in a newer solver I'm working on but one of the above tests did.

I can look at adding at least some of those with my up-coming PR for this new solver.