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 |