MAINT: fix most errors and warnings for `make html`

All that's left is some autodoc related warnings for scipy.odr
(from test_autodoc3.py)

2 files changed, 97 insertions, 2 deletions

diff --git a/examples/newton_krylov_preconditioning.py b/examples/newton_krylov_preconditioning.py
new file mode 100644
index 0000000..4736157
--- /dev/null
+++ b/examples/newton_krylov_preconditioning.py
@@ -0,0 +1,94 @@
+import numpy as np
+from scipy.optimize import root
+from scipy.sparse import spdiags, kron
+from scipy.sparse.linalg import spilu, LinearOperator
+from numpy import cosh, zeros_like, mgrid, zeros, eye
+
+# parameters
+nx, ny = 75, 75
+hx, hy = 1./(nx-1), 1./(ny-1)
+
+P_left, P_right = 0, 0
+P_top, P_bottom = 1, 0
+
+def get_preconditioner():
+    """Compute the preconditioner M"""
+    diags_x = zeros((3, nx))
+    diags_x[0,:] = 1/hx/hx
+    diags_x[1,:] = -2/hx/hx
+    diags_x[2,:] = 1/hx/hx
+    Lx = spdiags(diags_x, [-1,0,1], nx, nx)
+
+    diags_y = zeros((3, ny))
+    diags_y[0,:] = 1/hy/hy
+    diags_y[1,:] = -2/hy/hy
+    diags_y[2,:] = 1/hy/hy
+    Ly = spdiags(diags_y, [-1,0,1], ny, ny)
+
+    J1 = kron(Lx, eye(ny)) + kron(eye(nx), Ly)
+
+    # Now we have the matrix `J_1`. We need to find its inverse `M` --
+    # however, since an approximate inverse is enough, we can use
+    # the *incomplete LU* decomposition
+
+    J1_ilu = spilu(J1)
+
+    # This returns an object with a method .solve() that evaluates
+    # the corresponding matrix-vector product. We need to wrap it into
+    # a LinearOperator before it can be passed to the Krylov methods:
+
+    M = LinearOperator(shape=(nx*ny, nx*ny), matvec=J1_ilu.solve)
+    return M
+
+def solve(preconditioning=True):
+    """Compute the solution"""
+    count = [0]
+
+    def residual(P):
+        count[0] += 1
+
+        d2x = zeros_like(P)
+        d2y = zeros_like(P)
+
+        d2x[1:-1] = (P[2:]   - 2*P[1:-1] + P[:-2])/hx/hx
+        d2x[0]    = (P[1]    - 2*P[0]    + P_left)/hx/hx
+        d2x[-1]   = (P_right - 2*P[-1]   + P[-2])/hx/hx
+
+        d2y[:,1:-1] = (P[:,2:] - 2*P[:,1:-1] + P[:,:-2])/hy/hy
+        d2y[:,0]    = (P[:,1]  - 2*P[:,0]    + P_bottom)/hy/hy
+        d2y[:,-1]   = (P_top   - 2*P[:,-1]   + P[:,-2])/hy/hy
+
+        return d2x + d2y + 5*cosh(P).mean()**2
+
+    # preconditioner
+    if preconditioning:
+        M = get_preconditioner()
+    else:
+        M = None
+
+    # solve
+    guess = zeros((nx, ny), float)
+
+    sol = root(residual, guess, method='krylov',
+               options={'disp': True,
+                        'jac_options': {'inner_M': M}})
+    print('Residual', abs(residual(sol.x)).max())
+    print('Evaluations', count[0])
+
+    return sol.x
+
+def main():
+    sol = solve(preconditioning=True)
+
+    # visualize
+    import matplotlib.pyplot as plt
+    x, y = mgrid[0:1:(nx*1j), 0:1:(ny*1j)]
+    plt.clf()
+    plt.pcolor(x, y, sol)
+    plt.clim(0, 1)
+    plt.colorbar()
+    plt.show()
+
+
+if __name__ == "__main__":
+    main()
diff --git a/test_optimize.rst b/test_optimize.rst
index a09e0d7..0a55c28 100644
--- a/test_optimize.rst
+++ b/test_optimize.rst
@@ -29,7 +29,7 @@ The module contains:
 
 5. Multivariate equation system solvers (:func:`root`) using a variety of
    algorithms (e.g. hybrid Powell, Levenberg-Marquardt or large-scale
-   methods such as Newton-Krylov).
+   methods such as Newton-Krylov [KK]_).
 
 Below, several examples demonstrate their basic usage.
 
@@ -771,7 +771,8 @@ lot more depth to this topic than is shown here.
 
 .. rubric:: References
 
-Some further reading and related software:
+Some further reading and related software for solving large-scale problems
+( [PP]_, [AMG]_):
 
 .. [KK] D.A. Knoll and D.E. Keyes, "Jacobian-free Newton-Krylov methods",
         J. Comp. Phys. 193, 357 (2003).