"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "mpmath/matrices/eigen_symmetric.py" between
mpmath-0.19.tar.gz and mpmath-1.0.0.tar.gz

About: mpmath is a Python library for arbitrary-precision floating-point arithmetic.

eigen_symmetric.py  (mpmath-0.19):eigen_symmetric.py  (mpmath-1.0.0)
skipping to change at line 80 skipping to change at line 80
For a good introduction to Householder reflections, see also For a good introduction to Householder reflections, see also
Stoer, Bulirsch - Introduction to Numerical Analysis. Stoer, Bulirsch - Introduction to Numerical Analysis.
""" """
# note : the vector v of the i-th houshoulder reflector is stored in a[(i+1) :,i] # note : the vector v of the i-th houshoulder reflector is stored in a[(i+1) :,i]
# whereas v/<v,v> is stored in a[i,(i+1):] # whereas v/<v,v> is stored in a[i,(i+1):]
n = A.rows n = A.rows
for i in xrange(n - 1, 0, -1): for i in xrange(n - 1, 0, -1):
# scale the vector
# scale the vector
scale = 0 scale = 0
for k in xrange(0, i): for k in xrange(0, i):
scale += abs(A[k,i]) scale += abs(A[k,i])
scale_inv = 0 scale_inv = 0
if scale != 0: if scale != 0:
scale_inv = 1/scale scale_inv = 1/scale
# sadly there are floating point numbers not equal to zero whose reciprocal is infinity # sadly there are floating point numbers not equal to zero whose recipro cal is infinity
if i == 1 or scale == 0 or ctx.isinf(scale_inv): if i == 1 or scale == 0 or ctx.isinf(scale_inv):
E[i] = A[i-1,i] # nothing to do E[i] = A[i-1,i] # nothing to do
D[i] = 0 D[i] = 0
continue continue
# calculate parameters for housholder transfo rmation # calculate parameters for housholder transformation
H = 0 H = 0
for k in xrange(0, i): for k in xrange(0, i):
A[k,i] *= scale_inv A[k,i] *= scale_inv
H += A[k,i] * A[k,i] H += A[k,i] * A[k,i]
F = A[i-1,i] F = A[i-1,i]
G = ctx.sqrt(H) G = ctx.sqrt(H)
if F > 0: if F > 0:
G = -G G = -G
E[i] = scale * G E[i] = scale * G
H -= F * G H -= F * G
A[i-1,i] = F - G A[i-1,i] = F - G
F = 0 F = 0
# apply housholder transformation # apply housholder transformation
for j in xrange(0, i): for j in xrange(0, i):
if calc_ev: if calc_ev:
A[i,j] = A[j,i] / H A[i,j] = A[j,i] / H
G = 0 # calculate A*U G = 0 # calculate A*U
for k in xrange(0, j + 1): for k in xrange(0, j + 1):
G += A[k,j] * A[k,i] G += A[k,j] * A[k,i]
for k in xrange(j + 1, i): for k in xrange(j + 1, i):
G += A[j,k] * A[k,i] G += A[j,k] * A[k,i]
skipping to change at line 203 skipping to change at line 202
- Handbook for auto. comp., Vol II, Linear Algebra, p.212-226 (1971) - Handbook for auto. comp., Vol II, Linear Algebra, p.212-226 (1971)
For a good introduction to Householder reflections, see also For a good introduction to Householder reflections, see also
Stoer, Bulirsch - Introduction to Numerical Analysis. Stoer, Bulirsch - Introduction to Numerical Analysis.
""" """
n = A.rows n = A.rows
T[n-1] = 1 T[n-1] = 1
for i in xrange(n - 1, 0, -1): for i in xrange(n - 1, 0, -1):
# scale the vector # scale the vector
scale = 0 scale = 0
for k in xrange(0, i): for k in xrange(0, i):
scale += abs(ctx.re(A[k,i])) + abs(ctx.im(A[k,i])) scale += abs(ctx.re(A[k,i])) + abs(ctx.im(A[k,i]))
scale_inv = 0 scale_inv = 0
if scale != 0: if scale != 0:
scale_inv = 1 / scale scale_inv = 1 / scale
# sadly there are floating point number s not equal to zero whose reciprocal is infinity # sadly there are floating point numbers not equal to zero whose recipro cal is infinity
if scale == 0 or ctx.isinf(scale_inv): if scale == 0 or ctx.isinf(scale_inv):
E[i] = 0 E[i] = 0
D[i] = 0 D[i] = 0
T[i-1] = 1 T[i-1] = 1
continue continue
if i == 1: if i == 1:
F = A[i-1,i] F = A[i-1,i]
f = abs(F) f = abs(F)
E[i] = f E[i] = f
D[i] = 0 D[i] = 0
if f != 0: if f != 0:
T[i-1] = T[i] * F / f T[i-1] = T[i] * F / f
else: else:
T[i-1] = T[i] T[i-1] = T[i]
continue continue
# calculate parameters for housholder t ransformation # calculate parameters for housholder transformation
H = 0 H = 0
for k in xrange(0, i): for k in xrange(0, i):
A[k,i] *= scale_inv A[k,i] *= scale_inv
rr = ctx.re(A[k,i]) rr = ctx.re(A[k,i])
ii = ctx.im(A[k,i]) ii = ctx.im(A[k,i])
H += rr * rr + ii * ii H += rr * rr + ii * ii
F = A[i-1,i] F = A[i-1,i]
f = abs(F) f = abs(F)
skipping to change at line 255 skipping to change at line 254
E[i] = scale * G E[i] = scale * G
if f != 0: if f != 0:
F = F / f F = F / f
TZ = - T[i] * F # T[i-1]=-T[i]*F, but we need T[i-1] as temporary storage TZ = - T[i] * F # T[i-1]=-T[i]*F, but we need T[i-1] as temporary storage
G *= F G *= F
else: else:
TZ = -T[i] # T[i-1]=-T[i] TZ = -T[i] # T[i-1]=-T[i]
A[i-1,i] += G A[i-1,i] += G
F = 0 F = 0
# apply housholder transformation # apply housholder transformation
for j in xrange(0, i): for j in xrange(0, i):
A[i,j] = A[j,i] / H A[i,j] = A[j,i] / H
G = 0 # calculate A*U G = 0 # calculate A*U
for k in xrange(0, j + 1): for k in xrange(0, j + 1):
G += ctx.conj(A[k,j]) * A[k,i] G += ctx.conj(A[k,j]) * A[k,i]
for k in xrange(j + 1, i): for k in xrange(j + 1, i):
G += A[j,k] * A[k,i] G += A[j,k] * A[k,i]
skipping to change at line 278 skipping to change at line 277
HH = F / (2 * H) HH = F / (2 * H)
for j in xrange(0, i): # calculate reduced A for j in xrange(0, i): # calculate reduced A
F = A[j,i] F = A[j,i]
G = T[j] - HH * F # calculate Q G = T[j] - HH * F # calculate Q
T[j] = G T[j] = G
for k in xrange(0, j + 1): for k in xrange(0, j + 1):
A[k,j] -= ctx.conj(F) * T[k] + ctx.conj(G) * A[k,i] A[k,j] -= ctx.conj(F) * T[k] + ctx.conj(G) * A[k,i]
# as we use the lower left part for sto # as we use the lower left part for storage
rage # we have to use the transpose of the normal formula
# we have to use the transpose of the n
ormal formula
T[i-1] = TZ T[i-1] = TZ
D[i] = H D[i] = H
for i in xrange(1, n): # better for compatibility for i in xrange(1, n): # better for compatibility
E[i-1] = E[i] E[i-1] = E[i]
E[n-1] = 0 E[n-1] = 0
D[0] = 0 D[0] = 0
for i in xrange(0, n): for i in xrange(0, n):
skipping to change at line 866 skipping to change at line 865
elif isinstance(qtype, str): elif isinstance(qtype, str):
raise ValueError("unknown quadrature rule \"%s\"" % qtype) raise ValueError("unknown quadrature rule \"%s\"" % qtype)
elif not isinstance(qtype, str): elif not isinstance(qtype, str):
w = qtype(d, e) w = qtype(d, e)
else: else:
assert 0 assert 0
tridiag_eigen(ctx, d, e, z) tridiag_eigen(ctx, d, e, z)
for i in xrange(len(z)): for i in xrange(len(z)):
z[i] *= z[i] z[i] *= z[i]
z = z.transpose() z = z.transpose()
return (d, w * z) return (d, w * z)
################################################################################ ################## ################################################################################ ##################
################################################################################ ################## ################################################################################ ##################
################################################################################ ################## ################################################################################ ##################
def svd_r_raw(ctx, A, V = False, calc_u = False): def svd_r_raw(ctx, A, V = False, calc_u = False):
""" """
 End of changes. 10 change blocks. 
14 lines changed or deleted 11 lines changed or added

Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)