svdcmp.txt

VB常用数值算法集1 VB常用数值算法集1

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