unit Unit2;

interface
uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls,unit1, Forms, Dialogs;

procedure SVDCMP(A:matrx2; M, N:integer;var W:array of real; V:matrx2);
procedure SVBKSB(U:matrx2; W:array of real; V:matrx2;
                   M, N:integer; B:array of real;var X:array of real);
implementation
procedure SVDCMP(A:matrx2; M, N:integer;var W:array of real; V:matrx2);
label 1,2,3;
var
    RV1:array [0..100] of real;
    I,J,L,ITS,JJ,NM,K:integer;
    G,S,SCALE1,C,F,H,X,Y,Z,Sgn,ANORM,AAAAA:real;
begin
    If M < N Then
         ShowMessage('You must augment A with extra zero rows.');
    G:=0;
    SCALE1:=0;
    ANORM:=0;
    For I:=1 To N do
    begin
        L:=I + 1;
        RV1[I]:=SCALE1 * G;
        G:=0;
        S:=0;
        SCALE1:=0;
        If I <= M Then
        begin
            For K:=I To M do
                SCALE1:=SCALE1 + Abs(A[K, I]);
            If SCALE1 <> 0 Then
            begin
                For K:=I To M do
                begin
                    A[K, I]:=A[K, I] / SCALE1;
                    S:=S + A[K, I] * A[K, I];
                end;
                F:=A[I, I];
                If F > 0 Then
                    SGN:=1
                else
                    SGN:=-1;
                G:=-Sqrt(S) * Sgn;
                H:=F * G - S;
                A[I, I]:=F - G;
                If I <> N Then
                begin
                    For J:=L To N DO
                    begin
                        S:=0;
                        For K:=I To M do
                            S:=S + A[K, I] * A[K, J];
                        F:=S / H;
                        For K:=I To M do
                            A[K, J]:=A[K, J] + F * A[K, I];
                    end;
                end;
                For K:=I To M do
                    A[K, I]:=SCALE1 * A[K, I];
            end;
        end;
        W[I]:=SCALE1 * G;
        G:=0;
        S:=0;
        SCALE1:=0;
        If (I <= M) And (I <> N) Then
        begin
            For K:=L To N do
                SCALE1:=SCALE1 + Abs(A[I, K]);
            If SCALE1 <> 0 Then
            begin
                For K:=L To N do
                begin
                    A[I, K]:=A[I, K] / SCALE1;
                    S:=S + A[I, K] * A[I, K];
                end;
                F:=A[I, L];
                If F > 0 Then
                    SGN:=1
                else
                    SGN:=-1;
                G:=-Sqrt(S) * Sgn;
                H:=F * G - S;
                A[I, L]:=F - G;
                For K:=L To N do
                    RV1[K]:=A[I, K] / H;
                If I <> M Then
                begin
                    For J:=L To M do
                    begin
                        S:=0;
                        For K:=L To N do
                            S:=S + A[J, K] * A[I, K];
                        For K:=L To N do
                            A[J, K]:=A[J, K] + S * RV1[K];
                    end;
                end;
                For K:=L To N do
                    A[I, K]:=SCALE1 * A[I, K];
            end;
        end;
        If ANORM > Abs(W[I]) + Abs(RV1[I]) Then
            ANORM:=ANORM
        Else
            ANORM:=Abs(W[I]) + Abs(RV1[I]);
    end;
    For I:=N DownTo 1 do
    begin
        If I < N Then
        begin
            If G <> 0 Then
            begin
                For J:=L To N do
                    V[J, I]:=(A[I, J] / A[I, L]) / G;
                For J:=L To N do
                begin
                    S:=0;
                    For K:=L To N do
                        S:=S + A[I, K] * V[K, J];
                    For K:=L To N do
                        V[K, J]:=V[K, J] + S * V[K, I];
                end;
            end;
            For J:=L To N do
            begin
                V[I, J]:=0;
                V[J, I]:=0;
            end;
        end;
        V[I, I]:=1;
        G:=RV1[I];
        L:=I;
    end;
    For I:=N DownTo 1 do
    begin
        L:=I + 1;
        G:=W[I];
        If I < N Then
        begin
            For J:=L To N do
                A[I, J]:=0;
        end;
        If G <> 0 Then
        begin
            G:=1 / G;
            If I <> N Then
            begin
                For J:=L To N do
                begin
                    S:=0;
                    For K:=L To M do
                        S:=S + A[K, I] * A[K, J];
                    F:=(S / A[I, I]) * G;
                    For K:=I To M do
                        A[K, J]:=A[K, J] + F * A[K, I];
                end;
            end;
            For J:=I To M do
            begin
                A[J, I]:=A[J, I] * G;
            end;
        end
        Else
            For J:=I To M do
                A[J, I]:=0;
        A[I, I]:=A[I, I] + 1;
    end;
    For K:=N DownTo 1 do
    begin
        For ITS:=1 To 30 do
        begin
            For L:=K DownTo 1 do
            begin
                NM:=L - 1;
                If Abs(RV1[L]) + ANORM = ANORM Then GoTo 2;
                If Abs(W[NM]) + ANORM = ANORM Then GoTo 1;
            end;
1:          C:=0;
            S:=1;
            For I:=L To K do
            begin
                F:=S * RV1[I];
                If Abs(F) + ANORM <> ANORM Then
                begin
                    G:=W[I];
                    H:=Sqrt(F * F + G * G);
                    W[I]:=H;
                    H:=1 / H;
                    C:=(G * H);
                    S:=-(F * H);
                    For J:=1 To M do
                    begin
                        Y:=A[J, NM];
                        Z:=A[J, I];
                        A[J, NM]:=(Y * C) + (Z * S);
                        A[J, I]:=-(Y * S) + (Z * C);
                    end;
                end;
            end;
2:          Z:=W[K];
            If L = K Then
            begin
                If Z < 0 Then
                begin
                    W[K]:=-Z;
                    For J:=1 To N do
                        V[J, K]:=-V[J, K];
                end;
                GoTo 3;
            end;
            If ITS = 30 Then
                ShowMessage('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:=Sqrt(F * F + 1);
            If F > 0 Then
                    Sgn:=1
                else
                    Sgn:=-1;
            F:=((X - Z) * (X + Z) + H * ((Y / (F + Abs(G) * Sgn)) - H)) / X;
            C:=1;
            S:=1;
            For J:=L To NM do
            begin
                I:=J + 1;
                G:=RV1[I];
                Y:=W[I];
                H:=S * G;
                G:=G * C;
                Z:=Sqrt(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 do
                begin
                    X:=V[NM, J];
                    Z:=V[NM, I];
                    V[NM, J]:=(X * C) + (Z * S);
                    V[NM, I]:=-(X * S) + (Z * C);
                end;
                Z:=Sqrt(F * F + H * H);
                W[J]:=Z;
                If Z <> 0 Then
                begin
                    Z:=1 / Z;
                    C:=F * Z;
                    S:=H * Z;
                end;
                F:=(C * G) + (S * Y);
                X:=-(S * G) + (C * Y);
                For NM:=1 To M do
                begin
                    Y:=A[NM, J];
                    Z:=A[NM, I];
                    A[NM, J]:=(Y * C) + (Z * S);
                    A[NM, I]:=-(Y * S) + (Z * C);
                end;
            end;
            RV1[L]:=0;
            RV1[K]:=F;
            W[K]:=X;
        end;
3:  AAAAA:=1;
    end;
end;

procedure SVBKSB(U:matrx2; W:array of real; V:matrx2;
                   M, N:integer; B:array of real;var X:array of reaL);
var
    TMP:array[0..100] of real;
    I,J,JJ:integer;
    S:real;
begin
    For J:=1 To N do
    begin
        S:=0;
        If W[J] <> 0 Then
        begin
            For I:=1 To M do
                S:=S + U[I, J] * B[I];
            S:=S / W[J];
        end;
        TMP[J]:=S;
    end;
    For J:=1 To N do
    begin
        S:=0;
        For JJ:=1 To N do
            S:=S + V[J, JJ] * TMP[JJ];
        X[J]:=S;
    end;
end;
end.