Re: Help with translation: Fortran77 => Ada95

[tmoran@acm.org]

>Interesting that it's almost exactly twice as long as the Fortran.
  Would you prefer this version, shorter than the Fortran one?  I dropped
the comments that just showed the original Fortran code, tightened up
things that fit nicely on one line and didn't really need multiple lines,
pretty printed, and did a few other "personal preference" changes.

  Negative_Beta: exception;

  procedure Gauss(Alpha, Beta   : in     Vector;
                  Eps           : in     Real;
                  Zero, Weight  :    out Vector;
                  Ierr          :    out Integer) is
  --
  --         Input:  n - - the number of points in the Gaussian quadrature
  --                       formula; type integer
  --         Input:  alpha,beta - - arrays (of the same dimension) to be filled
  --                       with the values of  alpha(k-1), beta(k-1), k=1,2,
  --                       ...,n
  --                 eps - the relative accuracy desired in the nodes
  --                       and weights
  --
  --         Output: zero- array containing the Gaussian nodes (in increasing
  --                       order) zero(k)=x(k), k=1,2,...,n
  --                 weight - array containing the
  --                       Gaussian weights  weight(k)=w(k), k=1,2,...,n
  --                 ierr- an error flag equal to  0  on normal return,
  --                       equal to  i  if the QR algorithm does not
  --                       converge within 30 iterations on evaluating the
  --                       i-th eigenvalue.
  --         Exceptions:
  --                 Negative_Beta if one of the beta's is negative.
  --                 Constraint_Error if the input vectors don't match in
  --                        length, or have zero length.
  --
  --
    E : Vector(Alpha'range);   --  working space.
    J : Natural;
    I, K, M, Mml : Positive;
    P, G, R, S, C, F, B : Real;
  begin
    Ierr := 0;  -- hope for the best
    Zero(1) := Alpha(1);
    if Beta(1) < 0.0 then
      raise Negative_Beta;
    end if;
    Weight(1) := Beta(1);
    if Alpha'Length = 1 then
      return;
    end if;
    Weight(1) := 1.0;
    E(E'Last) := 0.0;
    for K in 2 .. Alpha'Last loop
      Zero(K) := Alpha(K);
      if Beta(K) < 0.0 then
        raise Negative_Beta;
      end if;
      E(K - 1) := Sqrt(Beta(K));
      Weight(K) := 0.0;
    end loop;
    for L in Alpha'range loop
      J := 0;
      -- Look for a small subdiagonal element.
      loop
          for Mm in L .. Alpha'Last loop
            M := Mm;
            exit when abs(E(M)) <= Eps * (abs(Zero(M)) + abs(Zero(M + 1)));
          end loop;
          P := Zero(L);
          exit when M = L;
          if J = 30 then
            Ierr := L;
            return;
          end if;
          J := J + 1;
          -- Form shift.
          G := (Zero(L + 1) - P) / (2.0 * E(L));
          R := Sqrt(G * G + 1.0);
          G := Zero(M) - P + E(L) / (G + Sign(R, G));
          S := 1.0;
          C := 1.0;
          P := 0.0;
          Mml := M - L;
          for I in reverse L .. M - 1 loop -- This hides outer I
            F := S * E(I);
            B := C * E(I);
            if abs(F) >= abs(G) then
              C := G / F;
              R := Sqrt(C * C + 1.0);
              E(I + 1) := F * R;
              S := 1.0 / R;
              C := C * S;
            else
              S := F / G;
              R := Sqrt(S * S + 1.0);
              E(I + 1) := G * R;
              C := 1.0 / R;
              S := S * C;
            end if;
            G := Zero(I + 1) - P;
            R := (Zero(I) - G) * S + 2.0 * C * B;
            P := S * R;
            Zero(I + 1) := G + P;
            G := C * R - B;
            -- Form first component of vector.
            F := Weight(I + 1);
            Weight(I + 1) := S * Weight(I) + C * F;
            Weight(I) := C * Weight(I) - S * F;
          end loop;
          Zero(L) := Zero(L) - P;
          E(L) := G;
          E(M) := 0.0;
      end loop;
    end loop;
    -- Order eigenvalues and eigenvectors.
    for Ii in 2 .. Alpha'Last loop
      I := Ii - 1;
      K := I;
      P := Zero(I);
      for J in Ii .. Alpha'Last loop -- This hides outer J
        if Zero(J) < P then
          K := J;
          P := Zero(J);
        end if;
      end loop;
      if K /= I then
        Zero(K) := Zero(I);
        Zero(I) := P;
        P := Weight(I);
        Weight(I) := Weight(K);
        Weight(K) := P;
      end if;
    end loop;
    for K in Weight'range loop -- This hides outer K
      Weight(K) := Beta(1) * Weight(K) * Weight(K);
    end loop;
  end Gauss;