[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
Added:
trunk/Lib/sandbox/multigrid/tests/
trunk/Lib/sandbox/multigrid/tests/test_relaxation.py
Modified:
trunk/Lib/sandbox/multigrid/relaxation.py
Log:
added relaxation unittests, additional relaxation methods.
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
+
+
Added: trunk/Lib/sandbox/multigrid/tests/test_relaxation.py
===================================================================
--- 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()
+
More information about the Scipy-svn
mailing list