?? svdcmp.txt
字號(hào):
Sub SVDCMP(A(), M, N, W(), V())
Dim RV1(100)
If M < N Then Print "You must augment A with extra zero rows."
G = 0#
SCALE1 = 0#
ANORM = 0#
For I = 1 To N
L = I + 1
RV1(I) = SCALE1 * G
G = 0#
S = 0#
SCALE1 = 0#
If I <= M Then
For K = I To M
SCALE1 = SCALE1 + Abs(A(K, I))
Next K
If SCALE1 <> 0# Then
For K = I To M
A(K, I) = A(K, I) / SCALE1
S = S + A(K, I) * A(K, I)
Next K
F = A(I, I)
G = -Sqr(S) * Sgn(F)
H = F * G - S
A(I, I) = F - G
If I <> N Then
For J = L To N
S = 0#
For K = I To M
S = S + A(K, I) * A(K, J)
Next K
F = S / H
For K = I To M
A(K, J) = A(K, J) + F * A(K, I)
Next K
Next J
End If
For K = I To M
A(K, I) = SCALE1 * A(K, I)
Next K
End If
End If
W(I) = SCALE1 * G
G = 0#
S = 0#
SCALE1 = 0#
If I <= M And I <> N Then
For K = L To N
SCALE1 = SCALE1 + Abs(A(I, K))
Next K
If SCALE1 <> 0# Then
For K = L To N
A(I, K) = A(I, K) / SCALE1
S = S + A(I, K) * A(I, K)
Next K
F = A(I, L)
G = -Sqr(S) * Sgn(F)
H = F * G - S
A(I, L) = F - G
For K = L To N
RV1(K) = A(I, K) / H
Next K
If I <> M Then
For J = L To M
S = 0#
For K = L To N
S = S + A(J, K) * A(I, K)
Next K
For K = L To N
A(J, K) = A(J, K) + S * RV1(K)
Next K
Next J
End If
For K = L To N
A(I, K) = SCALE1 * A(I, K)
Next K
End If
End If
If ANORM > Abs(W(I)) + Abs(RV1(I)) Then
ANORM = ANORM
Else
ANORM = Abs(W(I)) + Abs(RV1(I))
End If
Next I
For I = N To 1 Step -1
If I < N Then
If G <> 0# Then
For J = L To N
V(J, I) = (A(I, J) / A(I, L)) / G
Next J
For J = L To N
S = 0#
For K = L To N
S = S + A(I, K) * V(K, J)
Next K
For K = L To N
V(K, J) = V(K, J) + S * V(K, I)
Next K
Next J
End If
For J = L To N
V(I, J) = 0#
V(J, I) = 0#
Next J
End If
V(I, I) = 1#
G = RV1(I)
L = I
Next I
For I = N To 1 Step -1
L = I + 1
G = W(I)
If I < N Then
For J = L To N
A(I, J) = 0#
Next J
End If
If G <> 0# Then
G = 1# / G
If I <> N Then
For J = L To N
S = 0#
For K = L To M
S = S + A(K, I) * A(K, J)
Next K
F = (S / A(I, I)) * G
For K = I To M
A(K, J) = A(K, J) + F * A(K, I)
Next K
Next J
End If
For J = I To M
A(J, I) = A(J, I) * G
Next J
Else
For J = I To M
A(J, I) = 0#
Next J
End If
A(I, I) = A(I, I) + 1#
Next I
For K = N To 1 Step -1
For ITS = 1 To 30
For L = K To 1 Step -1
NM = L - 1
If Abs(RV1(L)) + ANORM = ANORM Then GoTo 2
If Abs(W(NM)) + ANORM = ANORM Then GoTo 1
Next L
1 C = 0#
S = 1#
For I = L To K
F = S * RV1(I)
If Abs(F) + ANORM <> ANORM Then
G = W(I)
H = Sqr(F * F + G * G)
W(I) = H
H = 1# / H
C = (G * H)
S = -(F * H)
For J = 1 To M
Y = A(J, NM)
Z = A(J, I)
A(J, NM) = (Y * C) + (Z * S)
A(J, I) = -(Y * S) + (Z * C)
Next J
End If
Next I
2 Z = W(K)
If L = K Then
If Z < 0# Then
W(K) = -Z
For J = 1 To N
V(J, K) = -V(J, K)
Next J
End If
GoTo 3
End If
If ITS = 30 Then Print "No convergence in 30 iterations"
X = W(L)
NM = K - 1
Y = W(NM)
G = RV1(NM)
H = RV1(K)
F = ((Y - Z) * (Y + Z) + (G - H) * (G + H)) / (2# * H * Y)
G = Sqr(F * F + 1#)
F = ((X - Z) * (X + Z) + H * ((Y / (F + Abs(G) * Sgn(F))) - H)) / X
C = 1#
S = 1#
For J = L To NM
I = J + 1
G = RV1(I)
Y = W(I)
H = S * G
G = G * C
Z = Sqr(F * F + H * H)
RV1(J) = Z
C = F / Z
S = H / Z
F = (X * C) + (G * S)
G = -(X * S) + (G * C)
H = Y * S
Y = Y * C
For NM = 1 To N
X = V(NM, J)
Z = V(NM, I)
V(NM, J) = (X * C) + (Z * S)
V(NM, I) = -(X * S) + (Z * C)
Next NM
Z = Sqr(F * F + H * H)
W(J) = Z
If Z <> 0# Then
Z = 1# / Z
C = F * Z
S = H * Z
End If
F = (C * G) + (S * Y)
X = -(S * G) + (C * Y)
For NM = 1 To M
Y = A(NM, J)
Z = A(NM, I)
A(NM, J) = (Y * C) + (Z * S)
A(NM, I) = -(Y * S) + (Z * C)
Next NM
Next J
RV1(L) = 0#
RV1(K) = F
W(K) = X
Next ITS
3 AAAAA = 1
Next K
End Sub?? 快捷鍵說(shuō)明
復(fù)制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號(hào)
Ctrl + =
減小字號(hào)
Ctrl + -