Overhaul the documentation

README.md

@@ -24,34 +24,13 @@ using NonlinearSolve, StaticArrays
 
 f(u, p) = u .* u .- 2
 u0 = @SVector[1.0, 1.0]
-probN = NonlinearProblem(f, u0)
-solver = solve(probN, NewtonRaphson(), abstol = 1e-9)
+prob = NonlinearProblem(f, u0)
+solver = solve(prob)
 
 ## Bracketing Methods
 
 f(u, p) = u .* u .- 2.0
 u0 = (1.0, 2.0) # brackets
-probB = IntervalNonlinearProblem(f, u0)
-sol = solve(probB, ITP())
+prob = IntervalNonlinearProblem(f, u0)
+sol = solve(prob)
 ```
-
-## v1.0 Breaking Release Highlights!
-
-v1.0 has been released for NonlinearSolve.jl, making it a decentralized solver library
-akin to DifferentialEquations.jl. For simple implementations of nonlinear solvers,
-you can now use SimpleNonlinearSolve.jl. `Falsi`, `Bisection`, and `NewtonRaphson`
-implementations designed for scalar and static vector inputs have all moved to the
-lower dependency version. NonlinearSolve.jl is thus designed for the larger scale
-more complex implementations, with `NewtonRaphson` now sporting support for
-LinearSolve.jl and soon SparseDiffTools.jl to allow for preconditioned Newton-Krylov and
-exploitation of sparsity. The two pieces will continue to grow in this direction,
-with NonlinearSolve.jl gaining more and more wrapped solver libraries and support
-for more complex methods, while SimpleNonlinearSolve.jl will keep a lower dependency
-version with implementations for small scale problems that do not need all of the
-extra tooling.
-
-Additionally, `NonlinearProblem` was split into `NonlinearProblem` and `IntervalNonlinearProblem`,
-i.e. the bracketing versions now have their own problem definition, rather than using
-a `Tuple` for `u0` in a `NonlinearProblem`. This helps for finding problem-algorithm
-pairing errors at type time and overall improves the documentation / makes the roles
-more clear.

docs/pages.jl

@@ -1,8 +1,12 @@
 # Put in a separate page so it can be used by SciMLDocs.jl
 
 pages = ["index.md",
-    "Tutorials" => Any["tutorials/nonlinear.md",
-        "tutorials/advanced.md",
+    "tutorials/getting_started.md",
+    "Tutorials" => Any[
+        "tutorials/code_optimization.md",
+        "tutorials/large_systems.md",
+        "tutorials/small_compile.md",
+        "tutorials/termination_conditions.md",
         "tutorials/iterator_interface.md"],
     "Basics" => Any["basics/NonlinearProblem.md",
         "basics/NonlinearFunctions.md",

docs/src/api/simplenonlinearsolve.md

@@ -23,7 +23,7 @@ These methods are suited for any general nonlinear root-finding problem , i.e. `
 ```@docs
 SimpleNewtonRaphson
 Broyden
-Halley
+SimpleHalley
 Klement
 SimpleTrustRegion
 SimpleDFSane

docs/src/basics/NonlinearProblem.md

@@ -1,4 +1,4 @@
-# Nonlinear Problems
+# [Nonlinear Problems](@id problems)
 
 ## The Three Types of Nonlinear Problems
 

docs/src/basics/NonlinearSolution.md

@@ -1,4 +1,4 @@
-# Nonlinear Solutions
+# [Nonlinear Solutions](@id solution)
 
 ```@docs
 SciMLBase.NonlinearSolution

docs/src/basics/TerminationCondition.md

@@ -1,4 +1,4 @@
-# Termination Conditions
+# [Termination Conditions](@id termination_condition)
 
 Provides a API to specify termination conditions for [`NonlinearProblem`](@ref) and
 [`SteadyStateProblem`](@ref). For details on the various termination modes, i.e.,

docs/src/solvers/BracketingSolvers.md

@@ -1,4 +1,4 @@
-# Interval Rootfinding Methods (Bracketing Solvers)
+# [Interval Rootfinding Methods (Bracketing Solvers)](@id bracketing)
 
 `solve(prob::IntervalNonlinearProblem,alg;kwargs)`
 

docs/src/solvers/NonlinearSystemSolvers.md

@@ -71,13 +71,19 @@ methods excel at small problems and problems defined with static arrays.
 
   - `SimpleNewtonRaphson()`: A simplified implementation of the Newton-Raphson method.
   - `Broyden()`: The classic Broyden's quasi-Newton method.
+  - `LBroyden()`: A low-memory Broyden implementation, similar to L-BFGS. This method is
+    common in machine learning contexts but is known to be unstable in comparison to many
+    other choices.
   - `Klement()`: A quasi-Newton method due to Klement. It's supposed to be more efficient
     than Broyden's method, and it seems to be in the cases that have been tried, but more
     benchmarking is required.
   - `SimpleTrustRegion()`: A dogleg trust-region Newton method. Improved globalizing stability
     for more robust fitting over basic Newton methods, though potentially with a cost.
   - `SimpleDFSane()`: A low-overhead implementation of the df-sane method for solving
     large-scale nonlinear systems of equations.
+  - `SimpleHalley()`: A low-overhead implementation of the Halley method. This is a higher order
+    method and thus can converge faster to low tolerances than a Newton method. Requires higher
+    order derivatives, so best used when automatic differentiation is available.
 
 !!! note
 
@@ -102,7 +108,6 @@ This is a wrapper package for importing solvers from NLsolve.jl into the SciML i
 
 Submethod choices for this algorithm include:
 
-  - `:fixedpoint`: Fixed-point iteration
   - `:anderson`: Anderson-accelerated fixed-point iteration
   - `:newton`: Classical Newton method with an optional line search
   - `:trust_region`: Trust region Newton method (the default choice)

docs/src/solvers/SteadyStateSolvers.md

@@ -1,4 +1,4 @@
-# Steady State Solvers
+# [Steady State Solvers](@id ss_solvers)
 
 `solve(prob::SteadyStateProblem,alg;kwargs)`
 

docs/src/tutorials/code_optimization.md

@@ -0,0 +1,58 @@
+# [Code Optimization for Solving Nonlinear Systems](@id code_optimization)
+
+## Code Optimization in Julia
+
+Before starting this tutorial, we recommend the reader to check out one of the
+many tutorials for optimization Julia code. The following is an incomplete
+list:
+
+  - [The Julia Performance Tips](https://docs.julialang.org/en/v1/manual/performance-tips/)
+  - [MIT 18.337 Course Notes on Optimizing Serial Code](https://mitmath.github.io/18337/lecture2/optimizing)
+  - [What scientists must know about hardware to write fast code](https://viralinstruction.com/posts/hardware/)
+
+User-side optimizations are important because, for sufficiently difficult problems,
+most time will be spent inside your `f` function, the function you are
+trying to solve. “Efficient solvers" are those that reduce the required
+number of `f` calls to hit the error tolerance. The main ideas for optimizing
+your nonlinear solver code, or any Julia function, are the following:
+
+  - Make it non-allocating
+  - Use StaticArrays for small arrays
+  - Use broadcast fusion
+  - Make it type-stable
+  - Reduce redundant calculations
+  - Make use of BLAS calls
+  - Optimize algorithm choice
+
+We'll discuss these strategies in the context of nonlinear solvers.
+Let's start with small systems.
+
+## Optimizing Nonlinear Solver Code for Small Systems
+
+```@example
+using NonlinearSolve, StaticArrays
+
+f(u, p) = u .* u .- p
+u0 = @SVector[1.0, 1.0]
+p = 2.0
+probN = NonlinearProblem(f, u0, p)
+sol = solve(probN, NewtonRaphson(), reltol = 1e-9)
+```
+
+## Using Jacobian Free Newton Krylov (JNFK) Methods
+
+If we want to solve the first example, without constructing the entire Jacobian
+
+```@example
+using NonlinearSolve, LinearSolve
+
+function f!(res, u, p)
+    @. res = u * u - p
+end
+u0 = [1.0, 1.0]
+p = 2.0
+prob = NonlinearProblem(f!, u0, p)
+
+linsolve = LinearSolve.KrylovJL_GMRES()
+sol = solve(prob, NewtonRaphson(; linsolve), reltol = 1e-9)
+```

docs/src/tutorials/getting_started.md

@@ -0,0 +1,148 @@
+# Getting Started with Nonlinear Rootfinding in Julia
+
+NonlinearSolve.jl is a system for solving rootfinding problems, i.e. finding
+the value $$u$$ such that $$f(u) = 0$$. In this tutorial we will go through
+the basics of NonlinearSolve.jl, demonstrating the core ideas and leading you
+to understanding the deeper parts of the documentation.
+
+## The Three Types of Nonlinear Systems
+
+There are three types of nonlinear systems:
+
+1. The "standard nonlinear system", i.e. the `NonlinearProblem`. This is a
+   system of equations with an initial condition where you want to satisfy
+   all equations simultaniously.
+2. The "interval rootfinding problem", i.e. the `IntervalNonlinearProblem`.
+   This is the case where you're given an interval `[a,b]` and need to find
+   where `f(u) = 0` for `u` inside the bounds.
+3. The "steady state problem", i.e. find the `u` such that `u' = f(u) = 0`.
+   While related to (1), it's not entirely the same because there's a uniquely
+   defined privledged root.
+4. The nonlinear least squares problem, which is an overconstrained nonlinear
+   system (i.e. more equations than states) which might not be satisfiable, i.e.
+   there may be no `u` such that `f(u) = 0`, and thus we find the `u` which
+   minimizes `||f(u)||` in the least squares sense.
+
+For now let's focus on the first two. The other two are covered in later tutorials,
+but from the first two we can show the general flow of the NonlinearSolve.jl package.
+
+## Problem Type 1: Solving Nonlinear Systems of Equations
+
+A nonlinear system $$f(u) = 0$$ is specified by defining a function `f(u,p)`,
+where `p` are the parameters of the system. For example, the following solves
+the vector equation $$f(u) = u^2 - p$$ for a vector of equations:
+
+```@example 1
+using NonlinearSolve
+
+f(u, p) = u .* u .- p
+u0 = [1.0, 1.0]
+p = 2.0
+prob = NonlinearProblem(f, u0, p)
+sol = solve(prob)
+```
+
+where `u0` is the initial condition for the rootfinder. Native NonlinearSolve.jl
+solvers use the given type of `u0` to determine the type used within the solver
+and the return. Note that the parameters `p` can be any type, but most are an
+AbstractArray for automatic differentiation.
+
+### Investigating the Solution
+
+To investigate the solution, one can look at the elements of the `NonlinearSolution`.
+The most important value is `sol.u`: this is the `u` that satisfies `f(u) = 0`. For example:
+
+```@example 1
+u = sol.u
+```
+
+```@example 1
+f(u, p)
+```
+
+This final value, the difference of the solution against zero, can also be found with `sol.resid`:
+
+```@example 1
+sol.resid
+```
+
+To know if the solution converged, or why the solution had not converged we can check the return
+code (`retcode`):
+
+```@example 1
+sol.retcode
+```
+
+There are multiple return codes which can mean the solve was successful, and thus we can use the
+general command `SciMLBase.successful_retcode` to check whether the solution process exited as
+intended:
+
+```@exmaple
+SciMLBase.successful_retcode(sol)
+```
+
+If we're curious about what it took to solve this equation, then you're in luck because all of the
+details can be found in `sol.stats`:
+
+```@example 1
+sol.stats
+```
+
+For more information on `NonlinearSolution`s, see the [`NonlinearSolution` manual page](@ref solution).
+
+### Interacting with the Solver Options
+
+While `sol = solve(prob)` worked for our case here, in many situations you may need to interact more
+deeply with how the solving process is done. First things first, you can specify the solver using the
+positional arguments. For example, let's set the solver to `TrustRegion()`:
+
+```@example 1
+solve(prob, TrustRegion())
+```
+
+For a complete list of solver choices, see [the nonlinear system solvers page](@ref nonlinearsystemsolvers).
+
+Next we can modify the tolerances. Here let's set some really low tolerances to force a tight solution:
+
+```@example 1
+solve(prob, TrustRegion(), reltol=1e-12, abstol=1e-12)
+```
+
+There are many more options for doing this configuring. Specifically for handling termination conditions,
+see the [Termination Conditions](@ref termination_condition) page for more details. And for more details on
+all of the available keyword arguments, see the [solver options](@ref solver_options) page.
+
+## Problem Type 2: Solving Interval Rootfinding Problems with Bracketing Methods
+
+For scalar rootfinding problems, bracketing methods exist in NonlinearSolve. The difference with bracketing
+methods is that with bracketing methods, instead of giving a `u0` initial condition, you pass a `uspan (a,b)`
+bracket in which the zero is expected to live. For example:
+
+```@example 1
+using NonlinearSolve
+f(u, p) = u * u - 2.0
+uspan = (1.0, 2.0) # brackets
+prob_int = IntervalNonlinearProblem(f, uspan)
+sol = solve(prob_int)
+```
+
+All of the same option handling from before works just as before, now just with different solver choices
+(see the [bracketing solvers](@ref bracketing) page for more details). For example, let's set the solver
+to `ITP()` and set a high absolute tolerance:
+
+```@example 1
+sol = solve(prob_int, ITP(), abstol = 0.01)
+```
+
+## Going Beyond the Basics: How to use the Documentation
+
+Congrats, you now know how to use the basics of NonlinearSolve.jl! However, there is so much more to
+see. Next check out:
+
+- [Some code optimization tricks to know about with NonlinearSolve.jl](@ref code_optimization)
+- [An iterator interface which lets you step through the solving process step by step](@ref iterator)
+- [How to handle large systems of equations efficiently](@ref large_systems)
+- [Ways to use NonlinearSolve.jl that is faster to startup and can statically compile](@ref fast_startup)
+- [More detailed termination conditions](@ref termination_conditions_tutorial)
+
+And also check out the rest of the manual.

docs/src/tutorials/iterator_interface.md

@@ -1,4 +1,8 @@
-# Nonlinear Solver Iterator Interface
+# [Nonlinear Solver Iterator Interface](@id iterator)
+
+!!! warn
+
+    This iterator interface will be expanded with a `step!` function soon!
 
 There is an iterator form of the nonlinear solver which mirrors the DiffEq integrator interface:
 

docs/src/tutorials/advanced.md → docs/src/tutorials/large_systems.md

@@ -1,4 +1,4 @@
-# Solving Large Ill-Conditioned Nonlinear Systems with NonlinearSolve.jl
+# [Solving Large Ill-Conditioned Nonlinear Systems with NonlinearSolve.jl](@id large_systems)
 
 This tutorial is for getting into the extra features of using NonlinearSolve.jl. Solving ill-conditioned nonlinear systems requires specializing the linear solver on properties of the Jacobian in order to cut down on the ``\mathcal{O}(n^3)`` linear solve and the ``\mathcal{O}(n^2)`` back-solves. This tutorial is designed to explain the advanced usage of NonlinearSolve.jl by solving the steady state stiff Brusselator partial differential equation (BRUSS) using NonlinearSolve.jl.