📄 d8r7.frm
字号:
VERSION 5.00
Begin VB.Form Form1
Caption = "Form1"
ClientHeight = 4245
ClientLeft = 60
ClientTop = 345
ClientWidth = 6075
LinkTopic = "Form1"
ScaleHeight = 4245
ScaleWidth = 6075
StartUpPosition = 3 'Windows Default
Begin VB.CommandButton Command1
Caption = "Command1"
Height = 375
Left = 3960
TabIndex = 0
Top = 3480
Width = 1215
End
End
Attribute VB_Name = "Form1"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
Private Sub Command1_Click()
'PROGRAM D8R7
'Driver for routine HQR
NP = 5
Dim A(5, 5), WR(5), WI(5)
A(1, 1) = 1#: A(1, 2) = 2#: A(1, 3) = 0#: A(1, 4) = 0#: A(1, 5) = 0#
A(2, 1) = -2#: A(2, 2) = 3#: A(2, 3) = 0#: A(2, 4) = 0#: A(2, 5) = 0#
A(3, 1) = 3#: A(3, 2) = 4#: A(3, 3) = 50#: A(3, 4) = 0#: A(3, 5) = 0#
A(4, 1) = -4#: A(4, 2) = 5#: A(4, 3) = -60#: A(4, 4) = 7#: A(4, 5) = 0#
A(5, 1) = -5#: A(5, 2) = 6#: A(5, 3) = -70#: A(5, 4) = 8#: A(5, 5) = -9#
Print Tab(5); " Matrix:"
Print Tab(5)
For I = 1 To NP
Print Tab(5);
For J = 1 To NP
Print Format$(A(I, J), "##.00"),
Next J
Print
Next I
Call BALANC(A(), NP)
Call ELMHES(A(), NP)
Call HQR(A(), NP, WR(), WI())
Print Tab(5)
Print Tab(5); "Eigenvalues:"
Print Tab(5)
Print Tab(5); " Real Imag."
For I = 1 To NP
Print Tab(5); Format$(WR(I), ".000000E+00");
Print Tab(19); Format$(WI(I), ".000000E+00")
Next I
End Sub
Sub HQR(A(), N, WR(), WI())
ANORM = Abs(A(1, 1))
For I = 2 To N
For J = I - 1 To N
ANORM = ANORM + Abs(A(I, J))
Next J
Next I
NN = N
T = 0#
1 If NN >= 1 Then
ITS = 0
2 For L = NN To 2 Step -1
S = Abs(A(L - 1, L - 1)) + Abs(A(L, L))
If S = 0# Then S = ANORM
If Abs(A(L, L - 1)) + S = S Then GoTo 3
Next L
L = 1
3 X = A(NN, NN)
If L = NN Then
WR(NN) = X + T
WI(NN) = 0#
NN = NN - 1
Else
Y = A(NN - 1, NN - 1)
W = A(NN, NN - 1) * A(NN - 1, NN)
If L = NN - 1 Then
P = 0.5 * (Y - X)
Q = P ^ 2 + W
Z = Sqr(Abs(Q))
X = X + T
If Q >= 0# Then
Z = P + Abs(Z) * Sgn(P)
WR(NN) = Z + X
WR(NN - 1) = WR(NN)
If Z <> 0# Then WR(NN) = X - W / Z
WI(NN) = 0#
WI(NN - 1) = 0#
Else
WR(NN) = X + P
WR(NN - 1) = WR(NN)
WI(NN) = Z
WI(NN - 1) = -Z
End If
NN = NN - 2
Else
If ITS = 30 Then Print " too many iterations "
If ITS = 10 Or ITS = 20 Then
T = T + X
For I = 1 To NN
A(I, I) = A(I, I) - X
Next I
S = Abs(A(NN, NN - 1)) + Abs(A(NN - 1, NN - 2))
X = 0.75 * S
Y = X
W = -0.4375 * S ^ 2
End If
ITS = ITS + 1
For M = NN - 2 To L Step -1
Z = A(M, M)
R = X - Z
S = Y - Z
P = (R * S - W) / A(M + 1, M) + A(M, M + 1)
Q = A(M + 1, M + 1) - Z - R - S
R = A(M + 2, M + 1)
S = Abs(P) + Abs(Q) + Abs(R)
P = P / S
Q = Q / S
R = R / S
If M = L Then GoTo 4
U = Abs(A(M, M - 1)) * (Abs(Q) + Abs(R))
BBB = Abs(A(M + 1, M + 1))
AAA = Abs(A(M - 1, M - 1)) + Abs(Z) + BBB
V = Abs(P) * AAA
If U + V = V Then GoTo 4
Next M
4 For I = M + 2 To NN
A(I, I - 2) = 0#
If I <> M + 2 Then A(I, I - 3) = 0#
Next I
For K = M To NN - 1
If K <> M Then
P = A(K, K - 1)
Q = A(K + 1, K - 1)
R = 0#
If K <> NN - 1 Then R = A(K + 2, K - 1)
X = Abs(P) + Abs(Q) + Abs(R)
If X <> 0# Then
P = P / X
Q = Q / X
R = R / X
End If
End If
S = Sqr(P ^ 2 + Q ^ 2 + R ^ 2) * Sgn(P)
If S <> 0# Then
If K = M Then
If L <> M Then A(K, K - 1) = -A(K, K - 1)
Else
A(K, K - 1) = -S * X
End If
P = P + S
X = P / S
Y = Q / S
Z = R / S
Q = Q / P
R = R / P
For J = K To NN
P = A(K, J) + Q * A(K + 1, J)
If K <> NN - 1 Then
P = P + R * A(K + 2, J)
A(K + 2, J) = A(K + 2, J) - P * Z
End If
A(K + 1, J) = A(K + 1, J) - P * Y
A(K, J) = A(K, J) - P * X
Next J
If NN > K + 3 Then
III = K + 3
Else
III = NN
End If
For I = L To III
P = X * A(I, K) + Y * A(I, K + 1)
If K <> NN - 1 Then
P = P + Z * A(I, K + 2)
A(I, K + 2) = A(I, K + 2) - P * R
End If
A(I, K + 1) = A(I, K + 1) - P * Q
A(I, K) = A(I, K) - P
Next I
End If
Next K
GoTo 2
End If
End If
GoTo 1
End If
End Sub
Sub BALANC(A(), N)
RADIX = 2#
SQRDX = RADIX ^ 2
RADIX = 2#
SQRDX = 4#
1 Last = 1
For I = 1 To N
C = 0#
R = 0#
For J = 1 To N
If J <> I Then
C = C + Abs(A(J, I))
R = R + Abs(A(I, J))
End If
Next J
If C <> 0# And R <> 0# Then
G = R / RADIX
F = 1#
S = C + R
2 If C < G Then
F = F * RADIX
C = C * SQRDX
GoTo 2
End If
G = R * RADIX
3 If C > G Then
F = F / RADIX
C = C / SQRDX
GoTo 3
End If
If (C + R) / F < 0.95 * S Then
Last = 0
G = 1# / F
For J = 1 To N
A(I, J) = A(I, J) * G
Next J
For J = 1 To N
A(J, I) = A(J, I) * F
Next J
End If
End If
Next I
If Last = 0 Then GoTo 1
End Sub
Sub ELMHES(A(), N)
If N > 2 Then
For M = 2 To N - 1
X = 0#
I = M
For J = M To N
If Abs(A(J, M - 1)) > Abs(X) Then
X = A(J, M - 1)
I = J
End If
Next J
If I <> M Then
For J = M - 1 To N
Y = A(I, J)
A(I, J) = A(M, J)
A(M, J) = Y
Next J
For J = 1 To N
Y = A(J, I)
A(J, I) = A(J, M)
A(J, M) = Y
Next J
End If
If X <> 0# Then
For I = M + 1 To N
Y = A(I, M - 1)
If Y <> 0# Then
Y = Y / X
A(I, M - 1) = Y
For J = M To N
A(I, J) = A(I, J) - Y * A(M, J)
Next J
For J = 1 To N
A(J, M) = A(J, M) + Y * A(J, I)
Next J
End If
Next I
End If
Next M
End If
End Sub
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -