1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198
|
#!/usr/bin/env python
#
""" Test functions for UMFPACK wrappers
"""
import warnings
import random
from numpy.testing import assert_array_almost_equal, dec, \
decorate_methods
from numpy.testing.utils import WarningManager
from scipy import rand, matrix, diag, eye
from scipy.sparse import csc_matrix, spdiags, SparseEfficiencyWarning
from scipy.sparse.linalg import linsolve
warnings.simplefilter('ignore',SparseEfficiencyWarning)
import numpy as np
try:
import scipy.sparse.linalg.dsolve.umfpack as um
except (ImportError, AttributeError):
_have_umfpack = False
else:
_have_umfpack = um.umfpack._um is not None
# Allow disabling of nose tests if umfpack not present
# See end of file for application
_umfpack_skip = dec.skipif(not _have_umfpack,
'UMFPACK appears not to be compiled')
class _DeprecationAccept:
def setUp(self):
self.mgr = WarningManager()
self.mgr.__enter__()
warnings.simplefilter("ignore", DeprecationWarning)
def tearDown(self):
self.mgr.__exit__()
class TestSolvers(_DeprecationAccept):
"""Tests inverting a sparse linear system"""
def test_solve_complex_without_umfpack(self):
"""Solve: single precision complex"""
linsolve.use_solver( useUmfpack = False )
a = self.a.astype('F')
b = self.b
x = linsolve.spsolve(a, b)
#print x
#print "Error: ", a*x-b
assert_array_almost_equal(a*x, b, decimal=4)
def test_solve_without_umfpack(self):
"""Solve: single precision"""
linsolve.use_solver( useUmfpack = False )
a = self.a.astype('f')
b = self.b
x = linsolve.spsolve(a, b.astype('f'))
#print x
#print "Error: ", a*x-b
assert_array_almost_equal(a*x, b, decimal=4)
def test_solve_complex_umfpack(self):
"""Solve with UMFPACK: double precision complex"""
linsolve.use_solver( useUmfpack = True )
a = self.a.astype('D')
b = self.b
x = linsolve.spsolve(a, b)
#print x
#print "Error: ", a*x-b
assert_array_almost_equal(a*x, b)
def test_solve_umfpack(self):
"""Solve with UMFPACK: double precision"""
linsolve.use_solver( useUmfpack = True )
a = self.a.astype('d')
b = self.b
x = linsolve.spsolve(a, b)
#print x
#print "Error: ", a*x-b
assert_array_almost_equal(a*x, b)
def test_solve_sparse_rhs(self):
"""Solve with UMFPACK: double precision, sparse rhs"""
linsolve.use_solver( useUmfpack = True )
a = self.a.astype('d')
b = csc_matrix( self.b )
x = linsolve.spsolve(a, b)
#print x
#print "Error: ", a*x-b
assert_array_almost_equal(a*x, self.b)
def test_factorized_umfpack(self):
"""Prefactorize (with UMFPACK) matrix for solving with multiple rhs"""
linsolve.use_solver( useUmfpack = True )
a = self.a.astype('d')
solve = linsolve.factorized( a )
x1 = solve( self.b )
assert_array_almost_equal(a*x1, self.b)
x2 = solve( self.b2 )
assert_array_almost_equal(a*x2, self.b2)
def test_factorized_without_umfpack(self):
"""Prefactorize matrix for solving with multiple rhs"""
linsolve.use_solver( useUmfpack = False )
a = self.a.astype('d')
solve = linsolve.factorized( a )
x1 = solve( self.b )
assert_array_almost_equal(a*x1, self.b)
x2 = solve( self.b2 )
assert_array_almost_equal(a*x2, self.b2)
def setUp(self):
self.a = spdiags([[1, 2, 3, 4, 5], [6, 5, 8, 9, 10]], [0, 1], 5, 5)
#print "The sparse matrix (constructed from diagonals):"
#print self.a
self.b = np.array([1, 2, 3, 4, 5])
self.b2 = np.array([5, 4, 3, 2, 1])
_DeprecationAccept.setUp(self)
class TestFactorization(_DeprecationAccept):
"""Tests factorizing a sparse linear system"""
def test_complex_lu(self):
"""Getting factors of complex matrix"""
umfpack = um.UmfpackContext("zi")
for A in self.complex_matrices:
umfpack.numeric(A)
(L,U,P,Q,R,do_recip) = umfpack.lu(A)
L = L.todense()
U = U.todense()
A = A.todense()
if not do_recip: R = 1.0/R
R = matrix(diag(R))
P = eye(A.shape[0])[P,:]
Q = eye(A.shape[1])[:,Q]
assert_array_almost_equal(P*R*A*Q,L*U)
def test_real_lu(self):
"""Getting factors of real matrix"""
umfpack = um.UmfpackContext("di")
for A in self.real_matrices:
umfpack.numeric(A)
(L,U,P,Q,R,do_recip) = umfpack.lu(A)
L = L.todense()
U = U.todense()
A = A.todense()
if not do_recip: R = 1.0/R
R = matrix(diag(R))
P = eye(A.shape[0])[P,:]
Q = eye(A.shape[1])[:,Q]
assert_array_almost_equal(P*R*A*Q,L*U)
def setUp(self):
random.seed(0) #make tests repeatable
self.real_matrices = []
self.real_matrices.append(spdiags([[1, 2, 3, 4, 5], [6, 5, 8, 9, 10]],
[0, 1], 5, 5) )
self.real_matrices.append(spdiags([[1, 2, 3, 4, 5], [6, 5, 8, 9, 10]],
[0, 1], 4, 5) )
self.real_matrices.append(spdiags([[1, 2, 3, 4, 5], [6, 5, 8, 9, 10]],
[0, 2], 5, 5) )
self.real_matrices.append(rand(3,3))
self.real_matrices.append(rand(5,4))
self.real_matrices.append(rand(4,5))
self.real_matrices = [csc_matrix(x).astype('d') for x \
in self.real_matrices]
self.complex_matrices = [x.astype(np.complex128)
for x in self.real_matrices]
_DeprecationAccept.setUp(self)
# Skip methods if umfpack not present
for cls in [TestSolvers, TestFactorization]:
decorate_methods(cls, _umfpack_skip)
if __name__ == "__main__":
import nose
nose.run(argv=['', __file__])
|