Merge pull request #11 from rgommers/donors

Remove "sponsored by Enthought" from banner

6 files changed, 99 insertions, 11 deletions

diff --git a/_theme/scipy/layout.html b/_theme/scipy/layout.html
index f0ab604..9d359d5 100644
--- a/_theme/scipy/layout.html
+++ b/_theme/scipy/layout.html
@@ -165,7 +165,7 @@
 <div class="container">
   <div class="top-scipy-org-logo-header">
     <a href="{{ pathto('index') }}">
-      <img style="border: 0;" alt="SciPy" src="{{ pathto('_static/img/scipy_org_logo.gif', 1) }}"></a>
+      <img style="border: 0;" alt="SciPy" src="{{ pathto('_static/img/scipy_org_logo.png', 1) }}"></a>
     </div>
   </div>
 </div>
diff --git a/_theme/scipy/static/img/scipy_org_logo.gif b/_theme/scipy/static/img/scipy_org_logo.gif
deleted file mode 100644
index ab2e3ac..0000000
--- a/_theme/scipy/static/img/scipy_org_logo.gif
+++ /dev/null
diff --git a/_theme/scipy/static/img/scipy_org_logo.png b/_theme/scipy/static/img/scipy_org_logo.png
new file mode 100644
index 0000000..a45cbe0
--- /dev/null
+++ b/_theme/scipy/static/img/scipy_org_logo.png
diff --git a/conf.py b/conf.py
index d8c50e0..4fbcabe 100644
--- a/conf.py
+++ b/conf.py
@@ -1,6 +1,6 @@
 needs_sphinx = '1.1'
 
-extensions = ['sphinx.ext.autodoc', 'sphinx.ext.pngmath', 'numpydoc',
+extensions = ['sphinx.ext.autodoc', 'sphinx.ext.imgmath', 'numpydoc',
               'sphinx.ext.intersphinx', 'sphinx.ext.coverage',
               'sphinx.ext.autosummary', 'matplotlib.sphinxext.plot_directive']
 
@@ -28,13 +28,6 @@ html_theme_options = {
                   ("http://docs.scipy.org/", "Docs")]
 }
 
-pngmath_latex_preamble = r"""
-\usepackage{color}
-\definecolor{textgray}{RGB}{51,51,51}
-\color{textgray}
-"""
-pngmath_use_preview = True
-pngmath_dvipng_args = ['-gamma 1.5', '-D 96', '-bg Transparent']
 
 #------------------------------------------------------------------------------
 # Plot style
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).