def qzdiv(stake, A, B, Q, Z):
'''
Takes U.T. matrices A, B, orthonormal matrices Q,Z, rearranges them
so that all cases of abs(B(i,i)/A(i,i))>stake are in lower right
corner, while preserving U.T. and orthonormal properties and Q'AZ' and Q'BZ'.
Parameters
----------
stake : number, dtype=float
A : array_like, dtype=float
An upper triangular matrix
B : array_like, dtype=float
An upper triangular matrix
Q : array_like, dtype=float
An orthonormal matrix from the QZ decomposition
Z : array_like, dtype=float
An orthonormal matrix from the QZ decomposition
Returns
-------
A : array_like, dtype=float
Rearranged A matrix
B : array_like, dtype=float
Rearranged B matrix
Q : array_like, dtype=float
Rearranged Q matrix
Z : array_like, dtype=float
Rearranged Z matrix
Notes
-----
Copyright: C.A. Sims, 1996, Yale University.
'''
n, jnk = A.shape
root = abs(vstack((np.diag(A), np.diag(B))).T)
tmp = (root[:,0]<1.e-13).astype(int)
root[:,0] = root[:,0]- tmp *(root[:,0]+root[:,1])
root[:,1] = root[:,1]/root[:,0]
for i in xrange(n,0,-1):
m=0
for j in xrange(i,0,-1):
if (root[j-1,1] > stake or root[j-1,1] < -.1):
m=j
break
if m==0:
print "qzdiv(): Inputs unchanged!"
return A, B, Q, Z
for k in xrange(m,i,1):
A, B, Q, Z = qzswitch(k,A,B,Q,Z)
tmp = root[k-1,1]
root[k-1,1] = root[k,1]
root[k,1] = tmp
return A, B, Q, Z
评论列表
文章目录