Solves linear systems with a tridiagonal coefficient matrix.
A solver that uses the Thomas algorithm (the Tridiagonal matrix
algorithm) for systems on the form
0
[b c ] [x] [d]
A.x = [a b c ] [x] = [d] = d.
[ a b c] [x] [d]
[ a b] [x] [d]
0
Parameters
----------
A : A sparse diagonal matrix (dia format) with shape n-by-p. The
coefficient matrix.
b : Numpy array, n-by-1. The right-hand-side vector.
Examples
--------
>>> import numpy as np
>>> import scipy.sparse as sparse
>>> import parsimony.utils.linalgs as linalgs
>>> np.random.seed(42)
>>>
>>> n = 10
>>> a = np.random.rand(n); a[-1] = 0.0
>>> b = np.random.rand(n)
>>> c = np.random.rand(n); c[0] = 0.0
>>> abc = np.vstack((a, b, c))
>>> A = sparse.dia_matrix((abc, [-1, 0, 1]), shape=(n, n))
>>> d = np.random.rand(n, 1)
>>>
>>> solver = linalgs.TridiagonalSolver()
>>> x = solver.solve(A, d)
>>> print x
[[ -1.84339326]
[ 4.62737333]
[-12.41571989]
[ 16.38029815]
[ 14.38143172]
[-14.58969243]
[ 6.21233944]
[ 1.34271395]
[ -1.63358708]
[ 4.88318651]]
>>> x_ = np.linalg.solve(A.toarray(), d)
>>> print x_
[[ -1.84339326]
[ 4.62737333]
[-12.41571989]
[ 16.38029815]
[ 14.38143172]
[-14.58969243]
[ 6.21233944]
[ 1.34271395]
[ -1.63358708]
[ 4.88318651]]
>>> np.linalg.norm(x - x_) < 5e-14
True
>>>
>>> import time
>>> n = 100
>>> a = np.random.rand(n); a[-1] = 0.0
>>> b = np.random.rand(n)
>>> c = np.random.rand(n); c[0] = 0.0
>>> abc = np.vstack((a, b, c))
>>> A = sparse.dia_matrix((abc, [-1, 0, 1]), shape=(n, n))
>>> d = np.random.rand(n, 1)
>>>
>>> t = time.time()
>>> x = solver.solve(A, d)
>>> print "Time:", time.time() - t # doctest: +SKIP
>>>
>>> t = time.time()
>>> x_ = np.linalg.solve(A.toarray(), d)
>>> print "Time:", time.time() - t # doctest: +SKIP
>>>
>>> np.linalg.norm(x - x_) < 5e-12
True
>>>
>>> n = 1000
>>> a = np.random.rand(n); a[-1] = 0.0
>>> b = np.random.rand(n)
>>> c = np.random.rand(n); c[0] = 0.0
>>> abc = np.vstack((a, b, c))
>>> A = sparse.dia_matrix((abc, [-1, 0, 1]), shape=(n, n))
>>> d = np.random.rand(n, 1)
>>>
>>> t = time.time()
>>> x = solver.solve(A, d)
>>> print "Time:", time.time() - t # doctest: +SKIP
>>>
>>> t = time.time()
>>> x_ = np.linalg.solve(A.toarray(), d)
>>> print "Time:", time.time() - t # doctest: +SKIP
>>>
>>> np.linalg.norm(x - x_) < 5e-9
True
- Overrides:
Solver.solve
|