# [Scipy-svn] r3161 - in trunk/Lib/sandbox/multigrid: . tests

scipy-svn@scip... scipy-svn@scip...
Thu Jul 12 03:08:18 CDT 2007

```Author: wnbell
Date: 2007-07-12 03:08:15 -0500 (Thu, 12 Jul 2007)
New Revision: 3161

trunk/Lib/sandbox/multigrid/tests/
trunk/Lib/sandbox/multigrid/tests/test_relaxation.py
Modified:
trunk/Lib/sandbox/multigrid/relaxation.py
Log:

Modified: trunk/Lib/sandbox/multigrid/relaxation.py
===================================================================
--- trunk/Lib/sandbox/multigrid/relaxation.py	2007-07-12 08:06:38 UTC (rev 3160)
+++ trunk/Lib/sandbox/multigrid/relaxation.py	2007-07-12 08:08:15 UTC (rev 3161)
@@ -1,7 +1,7 @@
import multigridtools
import numpy

-def gauss_seidel(A,x,b,iterations=1,sweep="forward"):
+def gauss_seidel(A,x,b,iterations=1,sweep='forward'):
"""
Perform Gauss-Seidel iteration on the linear system Ax=b

@@ -28,15 +28,19 @@

def jacobi(A,x,b,iterations=1,omega=1.0):
"""
-    Perform Gauss-Seidel iteration on the linear system Ax=b
+    Perform Jacobi iteration on the linear system Ax=b

-     Input:
+       x <- (1 - omega) x  +  omega * D^-1 (b - (A - D) x)
+
+    where D is the diagonal of A.
+
+    Input:
A - NxN csr_matrix
x - rank 1 ndarray of length N
b - rank 1 ndarray of length N
Optional:
iterations - number of iterations to perform (default: 1)
-         sweep      - slice of unknowns to relax (default: all in forward direction)
+         omega      - damping parameter (default: 1.0)
"""
sweep = slice(None)
(row_start,row_stop,row_step) = sweep.indices(A.shape[0])
@@ -54,3 +58,36 @@
omega)

+def polynomial_smoother(A,x,b,coeffs):
+    """
+    Apply a polynomial smoother to the system Ax=b
+
+    The smoother has the form:
+      x_new = x + p(A) (b - A*x)
+    where p(A) is a polynomial in A whose scalar coeffients
+    are specified (in decending order) by argument coeffs.
+
+    Eg.
+
+      Richardson iteration p(A) = c_0:
+         polynomial_smoother(A,x,b,[c_0])
+
+      Linear smoother p(A) = c_1*A + c_0:
+         polynomial_smoother(A,x,b,[c_1,c_0])
+
+      Quadratic smoother p(A) = c_2*A^2 + c_1*A + c_0:
+         polynomial_smoother(A,x,b,[c_2,c_1,c_0])
+
+
+    Note: Horner's Rule is applied to avoid computing A^k directly.
+    """
+
+    residual = (b - A*x)
+    h = coeffs[0]*residual
+
+    for c in coeffs[1:]:
+        h = c*residual + A*h
+
+    x += h
+
+

===================================================================
--- trunk/Lib/sandbox/multigrid/tests/test_relaxation.py	2007-07-12 08:06:38 UTC (rev 3160)
+++ trunk/Lib/sandbox/multigrid/tests/test_relaxation.py	2007-07-12 08:08:15 UTC (rev 3161)
@@ -0,0 +1,99 @@
+from numpy.testing import *
+
+import numpy
+import scipy
+from scipy import arange,ones,zeros,array,allclose
+from scipy.sparse import spdiags
+
+
+set_package_path()
+from scipy.multigrid.relaxation import polynomial_smoother,gauss_seidel,jacobi
+restore_path()
+
+
+class test_relaxation(NumpyTestCase):
+    def check_polynomial(self):
+        N  = 3
+        A  = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        x0 = arange(N).astype(numpy.float64)
+        x  = x0.copy()
+        b  = zeros(N)
+
+        r = (b - A*x0)
+        polynomial_smoother(A,x,b,[-1.0/3.0])
+
+        assert_almost_equal(x,x0-1.0/3.0*r)
+
+        x  = x0.copy()
+        polynomial_smoother(A,x,b,[0.2,-1])
+        assert_almost_equal(x,x0 + 0.2*A*r - r)
+
+        x  = x0.copy()
+        polynomial_smoother(A,x,b,[0.2,-1])
+        assert_almost_equal(x,x0 + 0.2*A*r - r)
+
+        x  = x0.copy()
+        polynomial_smoother(A,x,b,[-0.14285714,  1., -2.])
+        assert_almost_equal(x,x0 - 0.14285714*A*A*r + A*r - 2*r)
+
+
+    def check_gauss_seidel(self):
+        N = 1
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        x = arange(N).astype(numpy.float64)
+        b = zeros(N)
+        gauss_seidel(A,x,b)
+        assert_almost_equal(x,array([0]))
+
+        N = 3
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        x = arange(N).astype(numpy.float64)
+        b = zeros(N)
+        gauss_seidel(A,x,b)
+        assert_almost_equal(x,array([1.0/2.0,5.0/4.0,5.0/8.0]))
+
+        N = 1
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        x = arange(N).astype(numpy.float64)
+        b = zeros(N)
+        gauss_seidel(A,x,b,sweep='backward')
+        assert_almost_equal(x,array([0]))
+
+        N = 3
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        x = arange(N).astype(numpy.float64)
+        b = zeros(N)
+        gauss_seidel(A,x,b,sweep='backward')
+        assert_almost_equal(x,array([1.0/8.0,1.0/4.0,1.0/2.0]))
+
+        N = 1
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        x = arange(N).astype(numpy.float64)
+        b = array([10])
+        gauss_seidel(A,x,b)
+        assert_almost_equal(x,array([5]))
+
+        N = 3
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        x = arange(N).astype(numpy.float64)
+        b = array([10,20,30])
+        gauss_seidel(A,x,b)
+        assert_almost_equal(x,array([11.0/2.0,55.0/4,175.0/8.0]))
+
+
+        #forward and backward passes should give same result with x=ones(N),b=zeros(N)
+        N = 100
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        x = ones(N)
+        b = zeros(N)
+        gauss_seidel(A,x,b,iterations=200,sweep='forward')
+        resid1 = numpy.linalg.norm(A*x,2)
+        x = ones(N)
+        gauss_seidel(A,x,b,iterations=200,sweep='backward')
+        resid2 = numpy.linalg.norm(A*x,2)
+        self.assert_(resid1 < 0.01 and resid2 < 0.01)
+        self.assert_(allclose(resid1,resid2))
+
+if __name__ == '__main__':
+    NumpyTest().run()
+

```