for I := Lbound to Ubound do
      X[I] := B[I];
    for I := Lbound to Ubound do
      begin
        Ip := Index[I];
        Sum := X[Ip];
        X[Ip] := X[I];
        if K >= Lbound then
          for J := K to Pred(I) do
            { Sum := Sum - A[I,J] * X[J] }
            Sum := CSub(Sum, CMult(A[I,J], X[J]))
        else if CAbs(Sum) <> 0.0 then
          K := I;
        X[I] := Sum;
      end;
    for I := Ubound downto Lbound do
      begin
        Sum := X[I];
        if I < Ubound then
          for J := Succ(I) to Ubound do
            { Sum := Sum - A[I,J] * X[J]; }
            Sum := CSub(Sum, CMult(A[I,J], X[J]));
        { X[I] := Sum / A[I,I]; }
        X[I] := CDiv(Sum, A[I,I]);
      end;
  end;

  function SV_Decomp(A : TMatrix;
                     Lbound, Ubound1, Ubound2 : Integer;
                     S : TVector;
                     V : TMatrix) : Integer;
{ ----------------------------------------------------------------------
  This function is a translation of the EISPACK subroutine SVD

  This function determines the singular value decomposition
  A = U.S.V' of a real M by N rectangular matrix. Householder
  bidiagonalization and a variant of the QR algorithm are used.
  ----------------------------------------------------------------------
  This is a crude translation. Many of the original goto's
  have been kept!
  ---------------------------------------------------------------------- }

  var
    I, J, K, L, I1, K1, L1, Mn, Its : Integer;
    C, F, G, H, T, X, Y, Z, Tst1, Tst2, Scale : Float;
    R : TVector;

  label
    190, 210, 270, 290, 360, 390, 410, 430, 460,
    475, 490, 510, 520, 540, 565, 580, 650, 700;

  begin
    DimVector(R, Ubound2);
    Scale := 0.0;
    G := 0.0;
    X := 0.0;

    { Householder reduction to bidiagonal form }
    for I := Lbound to Ubound2 do
      begin
        L := I + 1;
        R[I] := Scale * G;
        G := 0.0;
        T := 0.0;
        Scale := 0.0;
        if I > Ubound1 then goto 210;

        for K := I to Ubound1 do
          Scale := Scale + Abs(A[K, I]);

        if Scale = 0.0 then goto 210;

        for K := I to Ubound1 do
          begin
            A[K, I] := A[K, I] / Scale;
            T := T + Sqr(A[K, I]);
          end;

        F := A[I, I];
        G := - DSgn(Sqrt(T), F);
        H := F * G - T;
        A[I, I] := F - G;
        if I = Ubound2 then goto 190;

        for J := L to Ubound2 do
          begin
            T := 0.0;
            for K := I to Ubound1 do
              T := T + A[K, I] * A[K, J];
            F := T / H;
            for K := I to Ubound1 do
              A[K, J] := A[K, J] + F * A[K, I];
          end;

190:    for K := I to Ubound1 do
          A[K, I] := Scale * A[K, I];

210:    S[I] := Scale * G;
        G := 0.0;
        T := 0.0;
        Scale := 0.0;
        if (I > Ubound1) or (I = Ubound2) then goto 290;

        for K := L to Ubound2 do
          Scale := Scale + Abs(A[I, K]);

        if Scale = 0.0 then goto 290;

        for K := L to Ubound2 do
          begin
            A[I, K] := A[I, K] / Scale;
            T := T + Sqr(A[I, K]);
          end;

        F := A[I, L];
        G := - DSgn(Sqrt(T), F);
        H := F * G - T;
        A[I, L] := F - G;

        for K := L to Ubound2 do
          R[K] := A[I, K] / H;

        if I = Ubound1 then goto 270;

        for J := L to Ubound1 do
          begin
            T := 0.0;
            for K := L to Ubound2 do
              T := T + A[J, K] * A[I, K];
            for K := L to Ubound2 do
              A[J, K] := A[J, K] + T * R[K];
          end;

270:    for K := L to Ubound2 do
          A[I, K] := Scale * A[I, K];

290:    X := Max(X, Abs(S[I]) + Abs(R[I]));
      end;

    { Accumulation of right-hand transformations }
    for I := Ubound2 downto Lbound do
      begin
        if I = Ubound2 then goto 390;
        if G = 0.0 then goto 360;

        for J := L to Ubound2 do
          { Double division avoids possible underflow }
          V[J, I] := (A[I, J] / A[I, L]) / G;

        for J := L to Ubound2 do
          begin
            T := 0.0;
            for K := L to Ubound2 do
              T := T + A[I, K] * V[K, J];
            for K := L to Ubound2 do
              V[K, J] := V[K, J] + T * V[K, I];
          end;

360:    for J := L to Ubound2 do
          begin
            V[I, J] := 0.0;
            V[J, I] := 0.0;
          end;

390:    V[I, I] := 1.0;
        G := R[I];
        L := I;
      end;


410:{ Accumulation of left-hand transformations }
    Mn := Min(Ubound1, Ubound2);

    for I := Mn downto Lbound do
      begin
        L := I + 1;
        G := S[I];
        if I = Ubound2 then goto 430;

        for J := L to Ubound2 do
          A[I, J] := 0.0;

430:    if G = 0.0 then goto 475;
        if I = Mn then goto 460;

        for J := L to Ubound2 do
          begin
            T := 0.0;

            for K := L to Ubound1 do
              T := T + A[K, I] * A[K, J];

            { Double division avoids possible underflow }
            F := (T / A[I, I]) / G;

            for K := I to Ubound1 do
              A[K, J] := A[K, J] + F * A[K, I];
          end;

460:    for J := I to Ubound1 do
          A[J, I] := A[J, I] / G;

        goto 490;

475:    for J := I to Ubound1 do
          A[J, I] := 0.0;

490:    A[I, I] := A[I, I] + 1.0;
      end;

510:{ Diagonalization of the bidiagonal form }
    Tst1 := X;
    for K := Ubound2 downto Lbound do
      begin
        K1 := K - 1;
        Its := 0;

520:    { Test for splitting }
        for L := K downto Lbound do
          begin
            L1 := L - 1;
            Tst2 := Tst1 + Abs(R[L]);
            if Tst2 = Tst1 then goto 565;
            { R[Lbound] is always zero, so there is no exit
                through the bottom of the loop  }
            Tst2 := Tst1 + Abs(S[L1]);
            if Tst2 = Tst1 then goto 540;
          end;

540:    { Cancellation of R[L] if L greater than 1 }
        C := 0.0;
        T := 1.0;

        for I := L to K do
          begin
            F := T * R[I];
            R[I] := C * R[I];
            Tst2 := Tst1 + Abs(F);
            if Tst2 = Tst1 then goto 565;
            G := S[I];
            H := Pythag(F, G);
            S[I] := H;
            C := G / H;
            T := - F / H;

            for J := Lbound to Ubound1 do
              begin
                Y := A[J, L1];
                Z := A[J, I];
                A[J, L1] := Y * C + Z * T;
                A[J, I] := - Y * T + Z * C;
              end;
           end;

565:    { Test for convergence }
        Z := S[K];
        if L = K then goto 650;

        if Its = 30 then
          begin
            SV_Decomp := - K;
            Exit;
          end;

        { Shift from bottom 2 by 2 minor }
        Its := Its + 1;
        X := S[L];
        Y := S[K1];
        G := R[K1];
        H := R[K];
        F := 0.5 * (((G + Z) / H) * ((G - Z) / Y) + Y / H - H / Y);
        G := Pythag(F, 1.0);
        F := X - (Z / X) * Z + (H / X) * (Y / (F + DSgn(G, F)) - H);

        { Next QR transformation }
        C := 1.0;
        T := 1.0;

        for I1 := L to K1 do
          begin
            I := I1 + 1;
            G := R[I];
            Y := S[I];
            H := T * G;
            G := C * G;
            Z := Pythag(F, H);
            R[I1] := Z;
            C := F / Z;
            T := H / Z;
            F := X * C + G * T;
            G := - X * T + G * C;
            H := Y * T;
            Y := Y * C;

            for J := Lbound to Ubound2 do
              begin
                X := V[J, I1];
                Z := V[J, I];
                V[J, I1] := X * C + Z * T;
                V[J, I] := - X * T + Z * C;
              end;

            Z := Pythag(F, H);
            S[I1] := Z;
            { Rotation can be arbitrary if Z is zero }
            if Z = 0.0 then goto 580;
            C := F / Z;
            T := H / Z;
580:        F := C * G + T * Y;
            X := - T * G + C * Y;

            for J := Lbound to Ubound1 do
              begin
                Y := A[J, I1];
                Z := A[J, I];
                A[J, I1] := Y * C + Z * T;
                A[J, I] := - Y * T + Z * C;
              end;
           end;

        R[L] := 0.0;
        R[K] := F;
        S[K] := X;
        goto 520;

650:    { Convergence }
        if Z >= 0.0 then goto 700;

        { S[K] is made non-negative }
        S[K] := - Z;

        for J := Lbound to Ubound2 do
          V[J, K] := - V[J, K];
700:  end;

    SV_Decomp := 0;
  end;

  procedure SV_SetZero(S : TVector; Lbound, Ubound : Integer; Tol : Float);
  var
    Threshold : Float;
    I         : Integer;
  begin
    Threshold := Tol * Max(S, Lbound, Ubound);
    for I := Lbound to Ubound do
      if S[I] < Threshold then S[I] := 0.0;
  end;

  procedure SV_Solve(U : TMatrix; S : TVector; V : TMatrix; B : TVector;
                     Lbound, Ubound1, Ubound2 : Integer;
                     X : TVector);
  var
    I, J, JJ : Integer;
    Sum      : Float;
    Tmp      : TVector;
  begin
    DimVector(Tmp, Ubound2);
    for J := Lbound to Ubound2 do
      begin
        Sum := 0.0;
        if S[J] > 0.0 then
          begin
            for I := Lbound to Ubound1 do
              Sum := Sum + U[I,J] * B[I];
            Sum := Sum / S[J];
          end;
        Tmp[J] := Sum;
      end;
    for J := Lbound to Ubound2 do
      begin
        Sum := 0.0;
        for JJ := Lbound to Ubound2 do
          Sum := Sum + V[J,JJ] * Tmp[JJ];
        X[J] := Sum;
      end;
  end;

  procedure SV_Approx(U : TMatrix; S : TVector; V : TMatrix;
                      Lbound, Ubound1, Ubound2 : Integer; A : TMatrix);
  var
    I, J, K : Integer;
  begin
    for I := Lbound to Ubound1 do
      for J := Lbound to Ubound2 do
        begin
          A[I,J] := 0.0;
          for K := Lbound to Ubound2 do
            if S[K] > 0.0 then
              A[I,J] := A[I,J] + U[I,K] * V[J,K];
        end;
  end;

  function QR_Decomp(A : TMatrix; Lbound, Ubound1, Ubound2 : Integer;
                     R : TMatrix) : Integer;
  var
    I, J,