Hey Everyone,
This patch implements linalg.solve_tridiagonal, whose omission was probably
a mistake.
Below is a patch to CVS pygsl/linalg.py, followed by a test program.
--- linalg.py.old 2006-04-17 16:10:53.000000000 -0700
+++ linalg.py 2006-04-17 16:14:06.000000000 -0700
@@ -729,6 +729,22 @@
# Tridiagonal Systems
#
+def solve_tridiag(diag, e, f, b):
+ """
+ returns x
+
+ This function solves the general N-by-N system A x = b where A is
+ tridiagonal. The form of A for the 4-by-4 case is shown below,
+
+ A = ( d_0 e_0 )
+ ( f_0 d_1 e_1 )
+ ( f_1 d_2 e_2 )
+ ( f_2 d_3 )
+ """
+ x = zeros(diag.shape, get_typecode(diag))
+ _gslwrap.gsl_linalg_solve_tridiag(diag, e, f, b, x)
+ return x
+
def solve_symm_tridiag(diag, e, b):
"""
returns x
#!/usr/bin/python
from __future__ import division
from pygsl import linalg
from pygsl._numobj import *
diag = ones((3,), Float)
subd = zeros((2,), Float)
supd = zeros((2,), Float)
rhs = array([1, 2, 3], Float)
trivial_1 = linalg.solve_symm_tridiag(diag, subd, rhs)
trivial_2 = linalg.solve_tridiag(diag, supd, subd, rhs)
print 'This vector should be zero:', (trivial_1 - trivial_2)
subd[1] = 100
asymm = linalg.solve_tridiag(diag, supd, subd, rhs)
print 'This vector shouldn\'t:', (trivial_1 - asymm)
--
Bruce J.A. Nourish <[EMAIL PROTECTED]> http://sps.la.asu.edu/~bjan/
>>> import this
-------------------------------------------------------
This SF.Net email is sponsored by xPML, a groundbreaking scripting language
that extends applications into web and mobile media. Attend the live webcast
and join the prime developer group breaking into this new coding territory!
http://sel.as-us.falkag.net/sel?cmd=lnk&kid=110944&bid=241720&dat=121642
_______________________________________________
pygsl-discuss mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/pygsl-discuss
