[Scipy-svn] r6321 - in trunk/scipy/linalg: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Fri Apr 9 22:54:35 EDT 2010
Author: warren.weckesser
Date: 2010-04-09 21:54:35 -0500 (Fri, 09 Apr 2010)
New Revision: 6321
Added:
trunk/scipy/linalg/tests/test_decomp_cholesky.py
Modified:
trunk/scipy/linalg/decomp_cholesky.py
trunk/scipy/linalg/tests/test_basic.py
trunk/scipy/linalg/tests/test_decomp.py
Log:
DOC+TEST: linalg: Add a docstring and tests for cho_solve_banded.
Modified: trunk/scipy/linalg/decomp_cholesky.py
===================================================================
--- trunk/scipy/linalg/decomp_cholesky.py 2010-04-10 01:55:23 UTC (rev 6320)
+++ trunk/scipy/linalg/decomp_cholesky.py 2010-04-10 02:54:35 UTC (rev 6321)
@@ -197,18 +197,39 @@
return c
-def cho_solve_banded((ab, lower), b, overwrite_b=False):
- """To be written..."""
+def cho_solve_banded((cb, lower), b, overwrite_b=False):
+ """Solve the linear equations A x = b, given the Cholesky factorization of A.
- ab = asarray_chkfinite(ab)
+ Parameters
+ ----------
+ (cb, lower) : tuple, (array, bool)
+ `cb` is the Cholesky factorization of A, as given by cholesky_banded.
+ `lower` must be the same value that was given to cholesky_banded.
+ b : array
+ Right-hand side
+ overwrite_b : bool
+ If True, the function will overwrite the values in `b`.
+
+ Returns
+ -------
+ x : array
+ The solution to the system A x = b
+
+ See also
+ --------
+ cholesky_banded : Cholesky factorization of a banded matrix
+
+ """
+
+ cb = asarray_chkfinite(cb)
b = asarray_chkfinite(b)
# Validate shapes.
- if ab.shape[-1] != b.shape[0]:
- raise ValueError("shapes of ab and b are not compatible.")
+ if cb.shape[-1] != b.shape[0]:
+ raise ValueError("shapes of cb and b are not compatible.")
- pbtrs, = get_lapack_funcs(('pbtrs',), (ab, b))
- x, info = pbtrs(ab, b, lower=lower, overwrite_b=overwrite_b)
+ pbtrs, = get_lapack_funcs(('pbtrs',), (cb, b))
+ x, info = pbtrs(cb, b, lower=lower, overwrite_b=overwrite_b)
if info > 0:
raise LinAlgError("%d-th leading minor not positive definite" % info)
if info < 0:
Modified: trunk/scipy/linalg/tests/test_basic.py
===================================================================
--- trunk/scipy/linalg/tests/test_basic.py 2010-04-10 01:55:23 UTC (rev 6320)
+++ trunk/scipy/linalg/tests/test_basic.py 2010-04-10 02:54:35 UTC (rev 6321)
@@ -335,69 +335,6 @@
assert_array_almost_equal(dot(a,x),b)
-class TestCholeskyBanded(TestCase):
-
- def test_upper_real(self):
- # Symmetric positive definite banded matrix `a`
- a = array([[4.0, 1.0, 0.0, 0.0],
- [1.0, 4.0, 0.5, 0.0],
- [0.0, 0.5, 4.0, 0.2],
- [0.0, 0.0, 0.2, 4.0]])
- # Banded storage form of `a`.
- ab = array([[-1.0, 1.0, 0.5, 0.2],
- [4.0, 4.0, 4.0, 4.0]])
- c = cholesky_banded(ab, lower=False)
- ufac = zeros_like(a)
- ufac[range(4),range(4)] = c[-1]
- ufac[(0,1,2),(1,2,3)] = c[0,1:]
- assert_array_almost_equal(a, dot(ufac.T, ufac))
-
- def test_upper_complex(self):
- # Hermitian positive definite banded matrix `a`
- a = array([[4.0, 1.0, 0.0, 0.0],
- [1.0, 4.0, 0.5, 0.0],
- [0.0, 0.5, 4.0, -0.2j],
- [0.0, 0.0, 0.2j, 4.0]])
- # Banded storage form of `a`.
- ab = array([[-1.0, 1.0, 0.5, -0.2j],
- [4.0, 4.0, 4.0, 4.0]])
- c = cholesky_banded(ab, lower=False)
- ufac = zeros_like(a)
- ufac[range(4),range(4)] = c[-1]
- ufac[(0,1,2),(1,2,3)] = c[0,1:]
- assert_array_almost_equal(a, dot(ufac.conj().T, ufac))
-
- def test_lower_real(self):
- # Symmetric positive definite banded matrix `a`
- a = array([[4.0, 1.0, 0.0, 0.0],
- [1.0, 4.0, 0.5, 0.0],
- [0.0, 0.5, 4.0, 0.2],
- [0.0, 0.0, 0.2, 4.0]])
- # Banded storage form of `a`.
- ab = array([[4.0, 4.0, 4.0, 4.0],
- [1.0, 0.5, 0.2, -1.0]])
- c = cholesky_banded(ab, lower=True)
- lfac = zeros_like(a)
- lfac[range(4),range(4)] = c[0]
- lfac[(1,2,3),(0,1,2)] = c[1,:3]
- assert_array_almost_equal(a, dot(lfac, lfac.T))
-
- def test_lower_complex(self):
- # Hermitian positive definite banded matrix `a`
- a = array([[4.0, 1.0, 0.0, 0.0],
- [1.0, 4.0, 0.5, 0.0],
- [0.0, 0.5, 4.0, -0.2j],
- [0.0, 0.0, 0.2j, 4.0]])
- # Banded storage form of `a`.
- ab = array([[4.0, 4.0, 4.0, 4.0],
- [1.0, 0.5, 0.2j, -1.0]])
- c = cholesky_banded(ab, lower=True)
- lfac = zeros_like(a)
- lfac[range(4),range(4)] = c[0]
- lfac[(1,2,3),(0,1,2)] = c[1,:3]
- assert_array_almost_equal(a, dot(lfac, lfac.conj().T))
-
-
class TestInv(TestCase):
def test_simple(self):
@@ -550,8 +487,6 @@
assert_array_almost_equal(x,direct_lstsq(a,b,1))
-
-
class TestPinv(TestCase):
def test_simple(self):
Modified: trunk/scipy/linalg/tests/test_decomp.py
===================================================================
--- trunk/scipy/linalg/tests/test_decomp.py 2010-04-10 01:55:23 UTC (rev 6320)
+++ trunk/scipy/linalg/tests/test_decomp.py 2010-04-10 02:54:35 UTC (rev 6321)
@@ -15,7 +15,8 @@
"""
import numpy as np
-from numpy.testing import *
+from numpy.testing import TestCase, assert_equal, assert_array_almost_equal, \
+ assert_array_equal, assert_raises, run_module_suite
from scipy.linalg import eig, eigvals, lu, svd, svdvals, cholesky, qr, \
schur, rsf2csf, lu_solve, lu_factor, solve, diagsvd, hessenberg, rq, \
@@ -790,55 +791,7 @@
def test_simple(self):
assert_array_almost_equal(diagsvd([1,0,0],3,3),[[1,0,0],[0,0,0],[0,0,0]])
-class TestCholesky(TestCase):
- def test_simple(self):
- a = [[8,2,3],[2,9,3],[3,3,6]]
- c = cholesky(a)
- assert_array_almost_equal(dot(transpose(c),c),a)
- c = transpose(c)
- a = dot(c,transpose(c))
- assert_array_almost_equal(cholesky(a,lower=1),c)
-
- def test_simple_complex(self):
- m = array([[3+1j,3+4j,5],[0,2+2j,2+7j],[0,0,7+4j]])
- a = dot(transpose(conjugate(m)),m)
- c = cholesky(a)
- a1 = dot(transpose(conjugate(c)),c)
- assert_array_almost_equal(a,a1)
- c = transpose(c)
- a = dot(c,transpose(conjugate(c)))
- assert_array_almost_equal(cholesky(a,lower=1),c)
-
- def test_random(self):
- n = 20
- for k in range(2):
- m = random([n,n])
- for i in range(n):
- m[i,i] = 20*(.1+m[i,i])
- a = dot(transpose(m),m)
- c = cholesky(a)
- a1 = dot(transpose(c),c)
- assert_array_almost_equal(a,a1)
- c = transpose(c)
- a = dot(c,transpose(c))
- assert_array_almost_equal(cholesky(a,lower=1),c)
-
- def test_random_complex(self):
- n = 20
- for k in range(2):
- m = random([n,n])+1j*random([n,n])
- for i in range(n):
- m[i,i] = 20*(.1+abs(m[i,i]))
- a = dot(transpose(conjugate(m)),m)
- c = cholesky(a)
- a1 = dot(transpose(conjugate(c)),c)
- assert_array_almost_equal(a,a1)
- c = transpose(c)
- a = dot(c,transpose(conjugate(c)))
- assert_array_almost_equal(cholesky(a,lower=1),c)
-
-
class TestQR(TestCase):
def test_simple(self):
Added: trunk/scipy/linalg/tests/test_decomp_cholesky.py
===================================================================
--- trunk/scipy/linalg/tests/test_decomp_cholesky.py (rev 0)
+++ trunk/scipy/linalg/tests/test_decomp_cholesky.py 2010-04-10 02:54:35 UTC (rev 6321)
@@ -0,0 +1,140 @@
+
+
+from numpy.testing import TestCase, assert_array_almost_equal
+
+from numpy import array, transpose, dot, conjugate, zeros_like
+from numpy.random import rand
+from scipy.linalg import cholesky, cholesky_banded, cho_solve_banded
+
+
+def random(size):
+ return rand(*size)
+
+
+class TestCholesky(TestCase):
+
+ def test_simple(self):
+ a = [[8,2,3],[2,9,3],[3,3,6]]
+ c = cholesky(a)
+ assert_array_almost_equal(dot(transpose(c),c),a)
+ c = transpose(c)
+ a = dot(c,transpose(c))
+ assert_array_almost_equal(cholesky(a,lower=1),c)
+
+ def test_simple_complex(self):
+ m = array([[3+1j,3+4j,5],[0,2+2j,2+7j],[0,0,7+4j]])
+ a = dot(transpose(conjugate(m)),m)
+ c = cholesky(a)
+ a1 = dot(transpose(conjugate(c)),c)
+ assert_array_almost_equal(a,a1)
+ c = transpose(c)
+ a = dot(c,transpose(conjugate(c)))
+ assert_array_almost_equal(cholesky(a,lower=1),c)
+
+ def test_random(self):
+ n = 20
+ for k in range(2):
+ m = random([n,n])
+ for i in range(n):
+ m[i,i] = 20*(.1+m[i,i])
+ a = dot(transpose(m),m)
+ c = cholesky(a)
+ a1 = dot(transpose(c),c)
+ assert_array_almost_equal(a,a1)
+ c = transpose(c)
+ a = dot(c,transpose(c))
+ assert_array_almost_equal(cholesky(a,lower=1),c)
+
+ def test_random_complex(self):
+ n = 20
+ for k in range(2):
+ m = random([n,n])+1j*random([n,n])
+ for i in range(n):
+ m[i,i] = 20*(.1+abs(m[i,i]))
+ a = dot(transpose(conjugate(m)),m)
+ c = cholesky(a)
+ a1 = dot(transpose(conjugate(c)),c)
+ assert_array_almost_equal(a,a1)
+ c = transpose(c)
+ a = dot(c,transpose(conjugate(c)))
+ assert_array_almost_equal(cholesky(a,lower=1),c)
+
+
+class TestCholeskyBanded(TestCase):
+ """Tests for cholesky_banded() and cho_solve_banded."""
+
+ def test_upper_real(self):
+ # Symmetric positive definite banded matrix `a`
+ a = array([[4.0, 1.0, 0.0, 0.0],
+ [1.0, 4.0, 0.5, 0.0],
+ [0.0, 0.5, 4.0, 0.2],
+ [0.0, 0.0, 0.2, 4.0]])
+ # Banded storage form of `a`.
+ ab = array([[-1.0, 1.0, 0.5, 0.2],
+ [4.0, 4.0, 4.0, 4.0]])
+ c = cholesky_banded(ab, lower=False)
+ ufac = zeros_like(a)
+ ufac[range(4),range(4)] = c[-1]
+ ufac[(0,1,2),(1,2,3)] = c[0,1:]
+ assert_array_almost_equal(a, dot(ufac.T, ufac))
+
+ b = array([0.0, 0.5, 4.2, 4.2])
+ x = cho_solve_banded((c, False), b)
+ assert_array_almost_equal(x, [0.0, 0.0, 1.0, 1.0])
+
+ def test_upper_complex(self):
+ # Hermitian positive definite banded matrix `a`
+ a = array([[4.0, 1.0, 0.0, 0.0],
+ [1.0, 4.0, 0.5, 0.0],
+ [0.0, 0.5, 4.0, -0.2j],
+ [0.0, 0.0, 0.2j, 4.0]])
+ # Banded storage form of `a`.
+ ab = array([[-1.0, 1.0, 0.5, -0.2j],
+ [4.0, 4.0, 4.0, 4.0]])
+ c = cholesky_banded(ab, lower=False)
+ ufac = zeros_like(a)
+ ufac[range(4),range(4)] = c[-1]
+ ufac[(0,1,2),(1,2,3)] = c[0,1:]
+ assert_array_almost_equal(a, dot(ufac.conj().T, ufac))
+
+ b = array([0.0, 0.5, 4.0-0.2j, 0.2j + 4.0])
+ x = cho_solve_banded((c, False), b)
+ assert_array_almost_equal(x, [0.0, 0.0, 1.0, 1.0])
+
+ def test_lower_real(self):
+ # Symmetric positive definite banded matrix `a`
+ a = array([[4.0, 1.0, 0.0, 0.0],
+ [1.0, 4.0, 0.5, 0.0],
+ [0.0, 0.5, 4.0, 0.2],
+ [0.0, 0.0, 0.2, 4.0]])
+ # Banded storage form of `a`.
+ ab = array([[4.0, 4.0, 4.0, 4.0],
+ [1.0, 0.5, 0.2, -1.0]])
+ c = cholesky_banded(ab, lower=True)
+ lfac = zeros_like(a)
+ lfac[range(4),range(4)] = c[0]
+ lfac[(1,2,3),(0,1,2)] = c[1,:3]
+ assert_array_almost_equal(a, dot(lfac, lfac.T))
+
+ b = array([0.0, 0.5, 4.2, 4.2])
+ x = cho_solve_banded((c, True), b)
+ assert_array_almost_equal(x, [0.0, 0.0, 1.0, 1.0])
+
+ def test_lower_complex(self):
+ # Hermitian positive definite banded matrix `a`
+ a = array([[4.0, 1.0, 0.0, 0.0],
+ [1.0, 4.0, 0.5, 0.0],
+ [0.0, 0.5, 4.0, -0.2j],
+ [0.0, 0.0, 0.2j, 4.0]])
+ # Banded storage form of `a`.
+ ab = array([[4.0, 4.0, 4.0, 4.0],
+ [1.0, 0.5, 0.2j, -1.0]])
+ c = cholesky_banded(ab, lower=True)
+ lfac = zeros_like(a)
+ lfac[range(4),range(4)] = c[0]
+ lfac[(1,2,3),(0,1,2)] = c[1,:3]
+ assert_array_almost_equal(a, dot(lfac, lfac.conj().T))
+
+ b = array([0.0, 0.5j, 3.8j, 3.8])
+ x = cho_solve_banded((c, True), b)
+ assert_array_almost_equal(x, [0.0, 0.0, 1.0j, 1.0])
More information about the Scipy-svn
mailing list