"Fossies" - the Fresh Open Source Software Archive  

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

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

eigen.py  (mpmath-1.0.0):eigen.py  (mpmath-1.1.0)
skipping to change at line 301 skipping to change at line 301
v = ctx.hypot(ctx.hypot(ctx.re(c), ctx.im(c)), ctx.hypot(ctx.re(s), ctx.im(s ))) v = ctx.hypot(ctx.hypot(ctx.re(c), ctx.im(c)), ctx.hypot(ctx.re(s), ctx.im(s )))
if v == 0: if v == 0:
v = 1 v = 1
c = 1 c = 1
s = 0 s = 0
else: else:
c /= v c /= v
s /= v s /= v
cc = ctx.conj(c)
cs = ctx.conj(s)
for k in xrange(n0, n): for k in xrange(n0, n):
# apply givens rotation from the left # apply givens rotation from the left
x = A[n0 ,k] x = A[n0 ,k]
y = A[n0+1,k] y = A[n0+1,k]
A[n0 ,k] = ctx.conj(c) * x + ctx.conj(s) * y A[n0 ,k] = cc * x + cs * y
A[n0+1,k] = -s * x + c * y A[n0+1,k] = c * y - s * x
for k in xrange(min(n1, n0+3)): for k in xrange(min(n1, n0+3)):
# apply givens rotation from the right # apply givens rotation from the right
x = A[k,n0 ] x = A[k,n0 ]
y = A[k,n0+1] y = A[k,n0+1]
A[k,n0 ] = c * x + s * y A[k,n0 ] = c * x + s * y
A[k,n0+1] = -ctx.conj(s) * x + ctx.conj(c) * y A[k,n0+1] = cc * y - cs * x
if not isinstance(Q, bool): if not isinstance(Q, bool):
for k in xrange(n): for k in xrange(n):
# eigenvectors # eigenvectors
x = Q[k,n0 ] x = Q[k,n0 ]
y = Q[k,n0+1] y = Q[k,n0+1]
Q[k,n0 ] = c * x + s * y Q[k,n0 ] = c * x + s * y
Q[k,n0+1] = -ctx.conj(s) * x + ctx.conj(c) * y Q[k,n0+1] = cc * y - cs * x
# chase the bulge # chase the bulge
for j in xrange(n0, n1 - 2): for j in xrange(n0, n1 - 2):
# calculate givens rotation # calculate givens rotation
c = A[j+1,j] c = A[j+1,j]
s = A[j+2,j] s = A[j+2,j]
v = ctx.hypot(ctx.hypot(ctx.re(c), ctx.im(c)), ctx.hypot(ctx.re(s), ctx. im(s))) v = ctx.hypot(ctx.hypot(ctx.re(c), ctx.im(c)), ctx.hypot(ctx.re(s), ctx. im(s)))
skipping to change at line 345 skipping to change at line 348
v = 1 v = 1
c = 1 c = 1
s = 0 s = 0
else: else:
A[j+1,j] = v A[j+1,j] = v
c /= v c /= v
s /= v s /= v
A[j+2,j] = 0 A[j+2,j] = 0
cc = ctx.conj(c)
cs = ctx.conj(s)
for k in xrange(j+1, n): for k in xrange(j+1, n):
# apply givens rotation from the left # apply givens rotation from the left
x = A[j+1,k] x = A[j+1,k]
y = A[j+2,k] y = A[j+2,k]
A[j+1,k] = ctx.conj(c) * x + ctx.conj(s) * y A[j+1,k] = cc * x + cs * y
A[j+2,k] = -s * x + c * y A[j+2,k] = c * y - s * x
for k in xrange(0, min(n1, j+4)): for k in xrange(0, min(n1, j+4)):
# apply givens rotation from the right # apply givens rotation from the right
x = A[k,j+1] x = A[k,j+1]
y = A[k,j+2] y = A[k,j+2]
A[k,j+1] = c * x + s * y A[k,j+1] = c * x + s * y
A[k,j+2] = -ctx.conj(s) * x + ctx.conj(c) * y A[k,j+2] = cc * y - cs * x
if not isinstance(Q, bool): if not isinstance(Q, bool):
for k in xrange(0, n): for k in xrange(0, n):
# eigenvectors # eigenvectors
x = Q[k,j+1] x = Q[k,j+1]
y = Q[k,j+2] y = Q[k,j+2]
Q[k,j+1] = c * x + s * y Q[k,j+1] = c * x + s * y
Q[k,j+2] = -ctx.conj(s) * x + ctx.conj(c) * y Q[k,j+2] = cc * y - cs * x
def hessenberg_qr(ctx, A, Q): def hessenberg_qr(ctx, A, Q):
""" """
This routine computes the Schur decomposition of an upper Hessenberg matrix A. This routine computes the Schur decomposition of an upper Hessenberg matrix A.
Given A, an unitary matrix Q is determined such that Given A, an unitary matrix Q is determined such that
Q' A Q = R and Q' Q = Q Q' = 1 Q' A Q = R and Q' Q = Q Q' = 1
where R is an upper right triangular matrix. Here ' denotes the hermitian where R is an upper right triangular matrix. Here ' denotes the hermitian
transpose (i.e. transposition and conjugation). transpose (i.e. transposition and conjugation).
 End of changes. 8 change blocks. 
12 lines changed or deleted 18 lines changed or added

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