Help with translation: Fortran77 => Ada95

Help with translation: Fortran77 => Ada95

[N&J]

Hi all,
Is anyone willing to help me translate a FORTRAN77 subroutine to Ada95? I
just can't understand one of those TERRIBLE and UNACCEPTABLE  Fortran "goto"

I attach the FORTRAN subroutine (It calculates the eigenvalues and the first
components of the eigenvectors of a symmetric tridiagonal matrix)

I can't cope with the loop starting at 105 and the "goto" to it before 240.

Thanks in advance,
John

-----SubRoutine gauss.f--------
c
c        Input:  n - - the number of points in the Gaussian quadrature
c                      formula; type integer
c                alpha,beta - - arrays of dimension  n  to be filled
c                      with the values of  alpha(k-1), beta(k-1), k=1,2,
c                      ...,n
c                eps - the relative accuracy desired in the nodes
c                      and weights
c
c        Output: zero- array of dimension  n  containing the Gaussian
c                      nodes (in increasing order)  zero(k)=x(k), k=1,2,
c                      ...,n
c                weight - array of dimension  n  containing the
c                      Gaussian weights  weight(k)=w(k), k=1,2,...,n
c                ierr- an error flag equal to  0  on normal return,
c                      equal to  i  if the QR algorithm does not
c                      converge within 30 iterations on evaluating the
c                      i-th eigenvalue, equal to  -1  if  n  is not in
c                      range, and equal to  -2  if one of the beta's is
c                      negative.
c
c The array  e  is needed for working space.
c
      dimension alpha(n),beta(n),zero(n),weight(n),e(n)
      if(n.lt.1) then
        ierr=-1
        return
      end if
      ierr=0
      zero(1)=alpha(1)
      if(beta(1).lt.0.) then
        ierr=-2
        return
      end if
      weight(1)=beta(1)
      if (n.eq.1) return
      weight(1)=1.
      e(n)=0.
      do 100 k=2,n
        zero(k)=alpha(k)
        if(beta(k).lt.0.) then
          ierr=-2
          return
        end if
        e(k-1)=sqrt(beta(k))
        weight(k)=0.
  100 continue
      do 240 l=1,n
        j=0
c
c Look for a small subdiagonal element.
c
  105   do 110 m=l,n
          if(m.eq.n) goto 120
          if(abs(e(m)).le.eps*(abs(zero(m))+abs(zero(m+1)))) goto 120
  110   continue
  120   p=zero(l)
        if(m.eq.l) goto 240
        if(j.eq.30) goto 400
        j=j+1
c
c Form shift.
c
        g=(zero(l+1)-p)/(2.*e(l))
        r=sqrt(g*g+1.)
        g=zero(m)-p+e(l)/(g+sign(r,g))
        s=1.
        c=1.
        p=0.
        mml=m-l
c
c For i=m-1 step -1 until l do ...
c
        do 200 ii=1,mml
          i=m-ii
          f=s*e(i)
          b=c*e(i)
          if(abs(f).lt.abs(g)) goto 150
          c=g/f
          r=sqrt(c*c+1.)
          e(i+1)=f*r
          s=1./r
          c=c*s
          goto 160
  150     s=f/g
          r=sqrt(s*s+1.)
          e(i+1)=g*r
          c=1./r
          s=s*c
  160     g=zero(i+1)-p
          r=(zero(i)-g)*s +2.*c*b
          p=s*r
          zero(i+1)=g+p
          g=c*r-b
c
c Form first component of vector.
c
          f=weight(i+1)
          weight(i+1)=s*weight(i)+c*f
          weight(i)=c*weight(i)-s*f
  200   continue
        zero(l)=zero(l)-p
        e(l)=g
        e(m)=0.
        goto 105
  240 continue
c
c Order eigenvalues and eigenvectors.
c
      do 300 ii=2,n
        i=ii-1
        k=i
        p=zero(i)
        do 260 j=ii,n
          if(zero(j).ge.p) goto 260
          k=j
          p=zero(j)
  260   continue
        if(k.eq.i) goto 300
        zero(k)=zero(i)
        zero(i)=p
        p=weight(i)
        weight(i)=weight(k)
        weight(k)=p
  300 continue
      do 310 k=1,n
        weight(k)=beta(1)*weight(k)*weight(k)
  310 continue
      return
c
c Set error - no convergence to an eigenvalue after 30 iterations.
c
  400 ierr=l
      return
      end

Re: Help with translation: Fortran77 => Ada95

[Warren W. Gay VE3WWG]

Gak... FORTRAN! Quick.. where is my wooden stake? ;-)

Warren.

N&J wrote:
> 
> Hi all,
> Is anyone willing to help me translate a FORTRAN77 subroutine to Ada95? I
> just can't understand one of those TERRIBLE and UNACCEPTABLE  Fortran "goto"
> 
> I attach the FORTRAN subroutine (It calculates the eigenvalues and the first
> components of the eigenvectors of a symmetric tridiagonal matrix)
> 
> I can't cope with the loop starting at 105 and the "goto" to it before 240.
> 
> Thanks in advance,
> John
> 
> -----SubRoutine gauss.f--------
> c
> c        Input:  n - - the number of points in the Gaussian quadrature
> c                      formula; type integer
> c                alpha,beta - - arrays of dimension  n  to be filled
> c                      with the values of  alpha(k-1), beta(k-1), k=1,2,
> c                      ...,n
> c                eps - the relative accuracy desired in the nodes
> c                      and weights
> c
> c        Output: zero- array of dimension  n  containing the Gaussian
> c                      nodes (in increasing order)  zero(k)=x(k), k=1,2,
> c                      ...,n
> c                weight - array of dimension  n  containing the
> c                      Gaussian weights  weight(k)=w(k), k=1,2,...,n
> c                ierr- an error flag equal to  0  on normal return,
> c                      equal to  i  if the QR algorithm does not
> c                      converge within 30 iterations on evaluating the
> c                      i-th eigenvalue, equal to  -1  if  n  is not in
> c                      range, and equal to  -2  if one of the beta's is
> c                      negative.
> c
> c The array  e  is needed for working space.
> c
>       dimension alpha(n),beta(n),zero(n),weight(n),e(n)
>       if(n.lt.1) then
>         ierr=-1
>         return
>       end if
>       ierr=0
>       zero(1)=alpha(1)
>       if(beta(1).lt.0.) then
>         ierr=-2
>         return
>       end if
>       weight(1)=beta(1)
>       if (n.eq.1) return
>       weight(1)=1.
>       e(n)=0.
>       do 100 k=2,n
>         zero(k)=alpha(k)
>         if(beta(k).lt.0.) then
>           ierr=-2
>           return
>         end if
>         e(k-1)=sqrt(beta(k))
>         weight(k)=0.
>   100 continue
>       do 240 l=1,n
>         j=0
> c
> c Look for a small subdiagonal element.
> c
>   105   do 110 m=l,n
>           if(m.eq.n) goto 120
>           if(abs(e(m)).le.eps*(abs(zero(m))+abs(zero(m+1)))) goto 120
>   110   continue
>   120   p=zero(l)
>         if(m.eq.l) goto 240
>         if(j.eq.30) goto 400
>         j=j+1
> c
> c Form shift.
> c
>         g=(zero(l+1)-p)/(2.*e(l))
>         r=sqrt(g*g+1.)
>         g=zero(m)-p+e(l)/(g+sign(r,g))
>         s=1.
>         c=1.
>         p=0.
>         mml=m-l
> c
> c For i=m-1 step -1 until l do ...
> c
>         do 200 ii=1,mml
>           i=m-ii
>           f=s*e(i)
>           b=c*e(i)
>           if(abs(f).lt.abs(g)) goto 150
>           c=g/f
>           r=sqrt(c*c+1.)
>           e(i+1)=f*r
>           s=1./r
>           c=c*s
>           goto 160
>   150     s=f/g
>           r=sqrt(s*s+1.)
>           e(i+1)=g*r
>           c=1./r
>           s=s*c
>   160     g=zero(i+1)-p
>           r=(zero(i)-g)*s +2.*c*b
>           p=s*r
>           zero(i+1)=g+p
>           g=c*r-b
> c
> c Form first component of vector.
> c
>           f=weight(i+1)
>           weight(i+1)=s*weight(i)+c*f
>           weight(i)=c*weight(i)-s*f
>   200   continue
>         zero(l)=zero(l)-p
>         e(l)=g
>         e(m)=0.
>         goto 105
>   240 continue
> c
> c Order eigenvalues and eigenvectors.
> c
>       do 300 ii=2,n
>         i=ii-1
>         k=i
>         p=zero(i)
>         do 260 j=ii,n
>           if(zero(j).ge.p) goto 260
>           k=j
>           p=zero(j)
>   260   continue
>         if(k.eq.i) goto 300
>         zero(k)=zero(i)
>         zero(i)=p
>         p=weight(i)
>         weight(i)=weight(k)
>         weight(k)=p
>   300 continue
>       do 310 k=1,n
>         weight(k)=beta(1)*weight(k)*weight(k)
>   310 continue
>       return
> c
> c Set error - no convergence to an eigenvalue after 30 iterations.
> c
>   400 ierr=l
>       return
>       end

-- 
Warren W. Gay VE3WWG
http://members.home.net/ve3wwg

Re: Help with translation: Fortran77 => Ada95

[Colin Paul Gloster]

Sorry, I am tired and am not going to look through all of that now nor
really think of it in terms of the mathematical problem domain which could
produce a more elegant solution...especially the latter since it may just be
neater to forget about porting the FORTRAN and just write a function
which does as required.

However I have a suggestion, bearing in mind abhorrance of goto. Depending
on how mistaken I am, in F77 the continue keyword on line 110  is part of the
looping construct of do on line 105. I.e. m starts at one and each time it is
incremented it is incremented by the value of n (so m is not necessarily just
increased by one each time) and the loop is to be run one hundred and ten times
with continue signalling each next iteration.

-- more FORTRAN 77
-- Look for a small subdiagonal element.
--
  105   do 110 m=l,n
          if(m.eq.n) goto 120
          if(abs(e(m)).le.eps*(abs(zero(m))+abs(zero(m+1)))) goto 120
  110   continue
  120   p=zero(l)
-- more FORTRAN 77

So in your Ada95 version you could try something like (tidy this up
in to proper compilable Ada when you have decided what you will go for)
-- more F77 to be ported
MAX_LOOP_ITERATIONS := 110;
loop_iterations_so_far := 1; -- or 0, you may change this in your translating
m := 1;
while loop_iterations_so_far <= MAX_LOOP_ITERATIONS  and  m /= n  and (
abs(e(m)) <= eps*(abs(zero(m))+abs(zero(m+1))) ) loop
               m := m + n;
end loop;
-- more F77 to be ported

Bye,
Colin Paul Gloster

Re: Help with translation: Fortran77 => Ada95

[Jeffrey Carter]

N&J wrote:
> 
> Hi all,
> Is anyone willing to help me translate a FORTRAN77 subroutine to Ada95? I
> just can't understand one of those TERRIBLE and UNACCEPTABLE  Fortran "goto"
> 
> I attach the FORTRAN subroutine (It calculates the eigenvalues and the first
> components of the eigenvectors of a symmetric tridiagonal matrix)
> 
> I can't cope with the loop starting at 105 and the "goto" to it before 240.
> 
> Thanks in advance,
> John
> 
> -----SubRoutine gauss.f--------
> c
> c        Input:  n - - the number of points in the Gaussian quadrature
> c                      formula; type integer
> c                alpha,beta - - arrays of dimension  n  to be filled
> c                      with the values of  alpha(k-1), beta(k-1), k=1,2,
> c                      ...,n
> c                eps - the relative accuracy desired in the nodes
> c                      and weights
> c
> c        Output: zero- array of dimension  n  containing the Gaussian
> c                      nodes (in increasing order)  zero(k)=x(k), k=1,2,
> c                      ...,n
> c                weight - array of dimension  n  containing the
> c                      Gaussian weights  weight(k)=w(k), k=1,2,...,n
> c                ierr- an error flag equal to  0  on normal return,
> c                      equal to  i  if the QR algorithm does not
> c                      converge within 30 iterations on evaluating the
> c                      i-th eigenvalue, equal to  -1  if  n  is not in
> c                      range, and equal to  -2  if one of the beta's is
> c                      negative.
> c
> c The array  e  is needed for working space.
> c
>       dimension alpha(n),beta(n),zero(n),weight(n),e(n)
>       if(n.lt.1) then
>         ierr=-1
>         return
>       end if
>       ierr=0
>       zero(1)=alpha(1)
>       if(beta(1).lt.0.) then
>         ierr=-2
>         return
>       end if
>       weight(1)=beta(1)
>       if (n.eq.1) return
>       weight(1)=1.
>       e(n)=0.
>       do 100 k=2,n
>         zero(k)=alpha(k)
>         if(beta(k).lt.0.) then
>           ierr=-2
>           return
>         end if
>         e(k-1)=sqrt(beta(k))
>         weight(k)=0.
>   100 continue
>       do 240 l=1,n
>         j=0
> c
> c Look for a small subdiagonal element.
> c
>   105   do 110 m=l,n
>           if(m.eq.n) goto 120
>           if(abs(e(m)).le.eps*(abs(zero(m))+abs(zero(m+1)))) goto 120
>   110   continue
>   120   p=zero(l)
>         if(m.eq.l) goto 240
>         if(j.eq.30) goto 400
>         j=j+1
> c
> c Form shift.
> c
>         g=(zero(l+1)-p)/(2.*e(l))
>         r=sqrt(g*g+1.)
>         g=zero(m)-p+e(l)/(g+sign(r,g))
>         s=1.
>         c=1.
>         p=0.
>         mml=m-l
> c
> c For i=m-1 step -1 until l do ...
> c
>         do 200 ii=1,mml
>           i=m-ii
>           f=s*e(i)
>           b=c*e(i)
>           if(abs(f).lt.abs(g)) goto 150
>           c=g/f
>           r=sqrt(c*c+1.)
>           e(i+1)=f*r
>           s=1./r
>           c=c*s
>           goto 160
>   150     s=f/g
>           r=sqrt(s*s+1.)
>           e(i+1)=g*r
>           c=1./r
>           s=s*c
>   160     g=zero(i+1)-p
>           r=(zero(i)-g)*s +2.*c*b
>           p=s*r
>           zero(i+1)=g+p
>           g=c*r-b
> c
> c Form first component of vector.
> c
>           f=weight(i+1)
>           weight(i+1)=s*weight(i)+c*f
>           weight(i)=c*weight(i)-s*f
>   200   continue
>         zero(l)=zero(l)-p
>         e(l)=g
>         e(m)=0.
>         goto 105
>   240 continue
> c
> c Order eigenvalues and eigenvectors.
> c
>       do 300 ii=2,n
>         i=ii-1
>         k=i
>         p=zero(i)
>         do 260 j=ii,n
>           if(zero(j).ge.p) goto 260
>           k=j
>           p=zero(j)
>   260   continue
>         if(k.eq.i) goto 300
>         zero(k)=zero(i)
>         zero(i)=p
>         p=weight(i)
>         weight(i)=weight(k)
>         weight(k)=p
>   300 continue
>       do 310 k=1,n
>         weight(k)=beta(1)*weight(k)*weight(k)
>   310 continue
>       return
> c
> c Set error - no convergence to an eigenvalue after 30 iterations.
> c
>   400 ierr=l
>       return
>       end

This at least compiles and has no Goto's (watch for wrapped lines):

with Ada.Numerics.Generic_Elementary_Functions;
procedure Make_It_Compile is
   type Real is digits 12; 
   type Vector is array (Positive range <>) of Real; 

   package Math_Stuff is new Ada.Numerics.Generic_Elementary_Functions
(Real);
   use Math_Stuff;

   Too_Short        : exception;
   Length_Mismatch  : exception;
   Negative_Beta    : exception;
   Did_Not_Converge : exception;

   procedure Gauss (Alpha  : in     Vector;
                    Beta   : in     Vector;
                    Eps    : in     Real;
                    Zero   :    out Vector;
                    Weight :    out Vector)
   is
      -- Translated from Fortran:
      --c
      --c        Input:  n - - the number of points in the Gaussian
quadrature
      --c                      formula; type integer
      --c                alpha,beta - - arrays of dimension  n  to be
filled
      --c                      with the values of  alpha(k-1),
beta(k-1), k=1,2,
      --c                      ...,n
      --c                eps - the relative accuracy desired in the
nodes
      --c                      and weights
      --c
      --c        Output: zero- array of dimension  n  containing the
Gaussian
      --c                      nodes (in increasing order) 
zero(k)=x(k), k=1,2,
      --c                      ...,n
      --c                weight - array of dimension  n  containing the
      --c                      Gaussian weights  weight(k)=w(k),
k=1,2,...,n
      --c                ierr- an error flag equal to  0  on normal
return,
      --c                      equal to  i  if the QR algorithm does not
      --c                      converge within 30 iterations on
evaluating the
      --c                      i-th eigenvalue, equal to  -1  if  n  is
not in
      --c                      range, and equal to  -2  if one of the
beta's is
      --c                      negative.
      --c
      --c The array  e  is needed for working space.
      --c

      -- Assume all arrays have lower bound of 1 (generalize later).
      
      -- Placeholder for Sign function
      function Sign (A : Real; B : Real) return Real is
         -- null;
      begin -- Sign
         return 1.0;
      end Sign;

      E   : Vector (Alpha'range);
      J   : Natural;
      M   : Positive;
      P   : Real;
      G   : Real;
      R   : Real;
      S   : Real;
      C   : Real;
      Mml : Positive;
      F   : Real;
      B   : Real;
      I   : Positive;
      K   : Positive;
   begin -- Gauss
      --      if(n.lt.1) then
      --        ierr=-1
      --        return
      --      end if
      if Alpha'Length < 1 then
         raise Too_Short;
      end if;

      if Alpha'Length /= Beta'Length or else
         Alpha'Length /= Zero'Length or else
         Alpha'Length /= Weight'Length
      then
         raise Length_Mismatch;
      end if;

      --      ierr=0

      Zero(1):=Alpha(1);

      --      if(beta(1).lt.0.) then
      --        ierr=-2
      --        return
      --      end if
      if Beta (1) < 0.0 then
         raise Negative_Beta;
      end if;

      Weight(1):=Beta(1);
      --      if (n.eq.1) return
      if Alpha'Length = 1 then
         return;
      end if;

      Weight(1):=1.0;
      E(E'Last):=0.0;
      --      do 100 k=2,n
      Do_100 : for K in 2 .. Alpha'Last loop
         Zero(K):=Alpha(K);
         --        if(beta(k).lt.0.) then
         --          ierr=-2
         --          return
         --        end if
         if Beta (K) < 0.0 then
            raise Negative_Beta;
         end if;

         E(K-1):=Sqrt(Beta(K));
         Weight(K):=0.0;
         --  100 continue
      end loop Do_100;

      --      do 240 l=1,n
      Do_240 : for L in Alpha'range loop
         J:=0;
         --c
         --c Look for a small subdiagonal element.
         --c
         Goto_105 : loop
            --  105   do 110 m=l,n
            Do_110 : for Mm in L .. Alpha'Last loop -- Determine M
               M := Mm;

               --          if(m.eq.n) goto 120
               exit Do_110 when M = Alpha'Last;

               --         
if(abs(e(m)).le.eps*(abs(zero(m))+abs(zero(m+1)))) goto 120
               exit Do_110 when abs E (M) <= Eps * (abs Zero (M) + abs
Zero (M + 1) );
               --  110   continue
            end loop Do_110;

            --  120   p=zero(l)
            P := Zero (L);

            --        if(m.eq.l) goto 240
            exit Goto_105 when M = L;

            --        if(j.eq.30) goto 400
            if J = 30 then
               raise Did_Not_Converge;
            end if;

            J:=J+1;
            --c
            --c Form shift.
            --c
            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;
            --c
            --c For i=m-1 step -1 until l do ...
            --c

            --        do 200 ii=1,mml
            Do_200 : for I in reverse L .. M - 1 loop -- This hides
outer I
               --          i=m-ii
               F:=S*E(I);
               B:=C*E(I);
               --          if(abs(f).lt.abs(g)) goto 150
               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;
                  --          goto 160
               else
                  --  150     s:=f/g
                  S := F / G;
                  R:=Sqrt(S*S+1.0);
                  E(I+1):=G*R;
                  C:=1.0/R;
                  S:=S*C;
               end if;

               --  160     g:=zero(i+1)-p
               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;
               --c
               --c Form first component of vector.
               --c
               F:=Weight(I+1);
               Weight(I+1):=S*Weight(I)+C*F;
               Weight(I):=C*Weight(I)-S*F;
               --  200   continue
            end loop Do_200;

            Zero(L):=Zero(L)-P;
            E(L):=G;
            E(M):=0.0;

            --        goto 105
         end loop Goto_105;
         --  240 continue
      end loop Do_240;
      --c
      --c Order eigenvalues and eigenvectors.
      --c
      --      do 300 ii=2,n
      Do_300 : for Ii in 2 .. Alpha'Last loop
         I:=Ii-1;
         K:=I;
         P:=Zero(I);

         --        do 260 j=ii,n
         Do_260 : for J in Ii .. Alpha'Last loop -- This hides outer J
            --          if(zero(j).ge.p) goto 260
            if Zero (J) < P then
               K:=J;
               P:=Zero(J);
            end if;
            --  260   continue
         end loop Do_260;

         --        if(k.eq.i) goto 300
         if K /= I then
            Zero(K):=Zero(I);
            Zero(I):=P;
            P:=Weight(I);
            Weight(I):=Weight(K);
            Weight(K):=P;
         end if;
         --  300 continue
      end loop Do_300;

      --      do 310 k=1,n
      Do_310 : for K in Weight'range loop -- This hides outer K
         Weight(K):=Beta(1)*Weight(K)*Weight(K);
         --  310 continue
      end loop Do_310;
      --      return
      --c
      --c Set error - no convergence to an eigenvalue after 30
iterations.
      --c
      --  400 ierr=l
      --      return
      --      end
   end Gauss;
begin -- Make_It_Compile;
   null;
end Make_It_Compile;

-- 
Jeff Carter
"You cheesy lot of second-hand electric donkey-bottom biters."
Monty Python & the Holy Grail

[SOLVED] Help with translation: Fortran77 => Ada95

[N&J]

Hi all,
I have been working on this translation for more than a day. Thanks to your
help, and from some personal e-mails I got and from a similar C ("goto"
free) program I found, I finally managed to translate it to Ada95 (I at
least think so!). Some minor modifications and testing remains. If anyone is
interested in the result please send me an e-mail and I will gladly send
him/her the code.

Thanks again,
John.

Re: Help with translation: Fortran77 => Ada95

[GianLuigi Piacentini]

    do label iteration_variable= starting value, ending value
[,                                                             
increment]
      statements to be looped
label last statement to be looped (usually continue)

so   do 100 k=2, n
     ...
100  continue

means loop statements in the ... place the number of times needed to k
to go from 2 to n (with increment 1)
each time k (the loop variable) is used in the loop, it takes its
current value (2, 3, .., n)

goto 120 is a way to jump out of the loop, control resuming at the
statement labelled 120.

goto's after an if seems quite clear, and in old code is used as a quick
and dirt way to skip some instructions (like old BASIC).

goto to the end of a loop is a way to jump to the end of a loop (usually
a continue) skipping some loop statement but remaining inside the loop.

Hoping this helps

G.L. Piacentini

Re: Help with translation: Fortran77 => Ada95

[Gary Scott]

Interesting that it's almost exactly twice as long as the Fortran.  Why
not just update the Fortran to Fortran 95 and make it a lot more clear
and a lot more comprehensible.


Jeffrey Carter wrote:
> 
> N&J wrote:
> >
> > Hi all,
> > Is anyone willing to help me translate a FORTRAN77 subroutine to Ada95? I
> > just can't understand one of those TERRIBLE and UNACCEPTABLE  Fortran "goto"
> >
> > I attach the FORTRAN subroutine (It calculates the eigenvalues and the first
> > components of the eigenvectors of a symmetric tridiagonal matrix)
> >
> > I can't cope with the loop starting at 105 and the "goto" to it before 240.
> >
> > Thanks in advance,
> > John
> >
> > -----SubRoutine gauss.f--------
> > c
> > c        Input:  n - - the number of points in the Gaussian quadrature
> > c                      formula; type integer
> > c                alpha,beta - - arrays of dimension  n  to be filled
> > c                      with the values of  alpha(k-1), beta(k-1), k=1,2,
> > c                      ...,n
> > c                eps - the relative accuracy desired in the nodes
> > c                      and weights
> > c
> > c        Output: zero- array of dimension  n  containing the Gaussian
> > c                      nodes (in increasing order)  zero(k)=x(k), k=1,2,
> > c                      ...,n
> > c                weight - array of dimension  n  containing the
> > c                      Gaussian weights  weight(k)=w(k), k=1,2,...,n
> > c                ierr- an error flag equal to  0  on normal return,
> > c                      equal to  i  if the QR algorithm does not
> > c                      converge within 30 iterations on evaluating the
> > c                      i-th eigenvalue, equal to  -1  if  n  is not in
> > c                      range, and equal to  -2  if one of the beta's is
> > c                      negative.
> > c
> > c The array  e  is needed for working space.
> > c
> >       dimension alpha(n),beta(n),zero(n),weight(n),e(n)
> >       if(n.lt.1) then
> >         ierr=-1
> >         return
> >       end if
> >       ierr=0
> >       zero(1)=alpha(1)
> >       if(beta(1).lt.0.) then
> >         ierr=-2
> >         return
> >       end if
> >       weight(1)=beta(1)
> >       if (n.eq.1) return
> >       weight(1)=1.
> >       e(n)=0.
> >       do 100 k=2,n
> >         zero(k)=alpha(k)
> >         if(beta(k).lt.0.) then
> >           ierr=-2
> >           return
> >         end if
> >         e(k-1)=sqrt(beta(k))
> >         weight(k)=0.
> >   100 continue
> >       do 240 l=1,n
> >         j=0
> > c
> > c Look for a small subdiagonal element.
> > c
> >   105   do 110 m=l,n
> >           if(m.eq.n) goto 120
> >           if(abs(e(m)).le.eps*(abs(zero(m))+abs(zero(m+1)))) goto 120
> >   110   continue
> >   120   p=zero(l)
> >         if(m.eq.l) goto 240
> >         if(j.eq.30) goto 400
> >         j=j+1
> > c
> > c Form shift.
> > c
> >         g=(zero(l+1)-p)/(2.*e(l))
> >         r=sqrt(g*g+1.)
> >         g=zero(m)-p+e(l)/(g+sign(r,g))
> >         s=1.
> >         c=1.
> >         p=0.
> >         mml=m-l
> > c
> > c For i=m-1 step -1 until l do ...
> > c
> >         do 200 ii=1,mml
> >           i=m-ii
> >           f=s*e(i)
> >           b=c*e(i)
> >           if(abs(f).lt.abs(g)) goto 150
> >           c=g/f
> >           r=sqrt(c*c+1.)
> >           e(i+1)=f*r
> >           s=1./r
> >           c=c*s
> >           goto 160
> >   150     s=f/g
> >           r=sqrt(s*s+1.)
> >           e(i+1)=g*r
> >           c=1./r
> >           s=s*c
> >   160     g=zero(i+1)-p
> >           r=(zero(i)-g)*s +2.*c*b
> >           p=s*r
> >           zero(i+1)=g+p
> >           g=c*r-b
> > c
> > c Form first component of vector.
> > c
> >           f=weight(i+1)
> >           weight(i+1)=s*weight(i)+c*f
> >           weight(i)=c*weight(i)-s*f
> >   200   continue
> >         zero(l)=zero(l)-p
> >         e(l)=g
> >         e(m)=0.
> >         goto 105
> >   240 continue
> > c
> > c Order eigenvalues and eigenvectors.
> > c
> >       do 300 ii=2,n
> >         i=ii-1
> >         k=i
> >         p=zero(i)
> >         do 260 j=ii,n
> >           if(zero(j).ge.p) goto 260
> >           k=j
> >           p=zero(j)
> >   260   continue
> >         if(k.eq.i) goto 300
> >         zero(k)=zero(i)
> >         zero(i)=p
> >         p=weight(i)
> >         weight(i)=weight(k)
> >         weight(k)=p
> >   300 continue
> >       do 310 k=1,n
> >         weight(k)=beta(1)*weight(k)*weight(k)
> >   310 continue
> >       return
> > c
> > c Set error - no convergence to an eigenvalue after 30 iterations.
> > c
> >   400 ierr=l
> >       return
> >       end
> 
> This at least compiles and has no Goto's (watch for wrapped lines):
> 
> with Ada.Numerics.Generic_Elementary_Functions;
> procedure Make_It_Compile is
>    type Real is digits 12;
>    type Vector is array (Positive range <>) of Real;
> 
>    package Math_Stuff is new Ada.Numerics.Generic_Elementary_Functions
> (Real);
>    use Math_Stuff;
> 
>    Too_Short        : exception;
>    Length_Mismatch  : exception;
>    Negative_Beta    : exception;
>    Did_Not_Converge : exception;
> 
>    procedure Gauss (Alpha  : in     Vector;
>                     Beta   : in     Vector;
>                     Eps    : in     Real;
>                     Zero   :    out Vector;
>                     Weight :    out Vector)
>    is
>       -- Translated from Fortran:
>       --c
>       --c        Input:  n - - the number of points in the Gaussian
> quadrature
>       --c                      formula; type integer
>       --c                alpha,beta - - arrays of dimension  n  to be
> filled
>       --c                      with the values of  alpha(k-1),
> beta(k-1), k=1,2,
>       --c                      ...,n
>       --c                eps - the relative accuracy desired in the
> nodes
>       --c                      and weights
>       --c
>       --c        Output: zero- array of dimension  n  containing the
> Gaussian
>       --c                      nodes (in increasing order)
> zero(k)=x(k), k=1,2,
>       --c                      ...,n
>       --c                weight - array of dimension  n  containing the
>       --c                      Gaussian weights  weight(k)=w(k),
> k=1,2,...,n
>       --c                ierr- an error flag equal to  0  on normal
> return,
>       --c                      equal to  i  if the QR algorithm does not
>       --c                      converge within 30 iterations on
> evaluating the
>       --c                      i-th eigenvalue, equal to  -1  if  n  is
> not in
>       --c                      range, and equal to  -2  if one of the
> beta's is
>       --c                      negative.
>       --c
>       --c The array  e  is needed for working space.
>       --c
> 
>       -- Assume all arrays have lower bound of 1 (generalize later).
> 
>       -- Placeholder for Sign function
>       function Sign (A : Real; B : Real) return Real is
>          -- null;
>       begin -- Sign
>          return 1.0;
>       end Sign;
> 
>       E   : Vector (Alpha'range);
>       J   : Natural;
>       M   : Positive;
>       P   : Real;
>       G   : Real;
>       R   : Real;
>       S   : Real;
>       C   : Real;
>       Mml : Positive;
>       F   : Real;
>       B   : Real;
>       I   : Positive;
>       K   : Positive;
>    begin -- Gauss
>       --      if(n.lt.1) then
>       --        ierr=-1
>       --        return
>       --      end if
>       if Alpha'Length < 1 then
>          raise Too_Short;
>       end if;
> 
>       if Alpha'Length /= Beta'Length or else
>          Alpha'Length /= Zero'Length or else
>          Alpha'Length /= Weight'Length
>       then
>          raise Length_Mismatch;
>       end if;
> 
>       --      ierr=0
> 
>       Zero(1):=Alpha(1);
> 
>       --      if(beta(1).lt.0.) then
>       --        ierr=-2
>       --        return
>       --      end if
>       if Beta (1) < 0.0 then
>          raise Negative_Beta;
>       end if;
> 
>       Weight(1):=Beta(1);
>       --      if (n.eq.1) return
>       if Alpha'Length = 1 then
>          return;
>       end if;
> 
>       Weight(1):=1.0;
>       E(E'Last):=0.0;
>       --      do 100 k=2,n
>       Do_100 : for K in 2 .. Alpha'Last loop
>          Zero(K):=Alpha(K);
>          --        if(beta(k).lt.0.) then
>          --          ierr=-2
>          --          return
>          --        end if
>          if Beta (K) < 0.0 then
>             raise Negative_Beta;
>          end if;
> 
>          E(K-1):=Sqrt(Beta(K));
>          Weight(K):=0.0;
>          --  100 continue
>       end loop Do_100;
> 
>       --      do 240 l=1,n
>       Do_240 : for L in Alpha'range loop
>          J:=0;
>          --c
>          --c Look for a small subdiagonal element.
>          --c
>          Goto_105 : loop
>             --  105   do 110 m=l,n
>             Do_110 : for Mm in L .. Alpha'Last loop -- Determine M
>                M := Mm;
> 
>                --          if(m.eq.n) goto 120
>                exit Do_110 when M = Alpha'Last;
> 
>                --
> if(abs(e(m)).le.eps*(abs(zero(m))+abs(zero(m+1)))) goto 120
>                exit Do_110 when abs E (M) <= Eps * (abs Zero (M) + abs
> Zero (M + 1) );
>                --  110   continue
>             end loop Do_110;
> 
>             --  120   p=zero(l)
>             P := Zero (L);
> 
>             --        if(m.eq.l) goto 240
>             exit Goto_105 when M = L;
> 
>             --        if(j.eq.30) goto 400
>             if J = 30 then
>                raise Did_Not_Converge;
>             end if;
> 
>             J:=J+1;
>             --c
>             --c Form shift.
>             --c
>             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;
>             --c
>             --c For i=m-1 step -1 until l do ...
>             --c
> 
>             --        do 200 ii=1,mml
>             Do_200 : for I in reverse L .. M - 1 loop -- This hides
> outer I
>                --          i=m-ii
>                F:=S*E(I);
>                B:=C*E(I);
>                --          if(abs(f).lt.abs(g)) goto 150
>                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;
>                   --          goto 160
>                else
>                   --  150     s:=f/g
>                   S := F / G;
>                   R:=Sqrt(S*S+1.0);
>                   E(I+1):=G*R;
>                   C:=1.0/R;
>                   S:=S*C;
>                end if;
> 
>                --  160     g:=zero(i+1)-p
>                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;
>                --c
>                --c Form first component of vector.
>                --c
>                F:=Weight(I+1);
>                Weight(I+1):=S*Weight(I)+C*F;
>                Weight(I):=C*Weight(I)-S*F;
>                --  200   continue
>             end loop Do_200;
> 
>             Zero(L):=Zero(L)-P;
>             E(L):=G;
>             E(M):=0.0;
> 
>             --        goto 105
>          end loop Goto_105;
>          --  240 continue
>       end loop Do_240;
>       --c
>       --c Order eigenvalues and eigenvectors.
>       --c
>       --      do 300 ii=2,n
>       Do_300 : for Ii in 2 .. Alpha'Last loop
>          I:=Ii-1;
>          K:=I;
>          P:=Zero(I);
> 
>          --        do 260 j=ii,n
>          Do_260 : for J in Ii .. Alpha'Last loop -- This hides outer J
>             --          if(zero(j).ge.p) goto 260
>             if Zero (J) < P then
>                K:=J;
>                P:=Zero(J);
>             end if;
>             --  260   continue
>          end loop Do_260;
> 
>          --        if(k.eq.i) goto 300
>          if K /= I then
>             Zero(K):=Zero(I);
>             Zero(I):=P;
>             P:=Weight(I);
>             Weight(I):=Weight(K);
>             Weight(K):=P;
>          end if;
>          --  300 continue
>       end loop Do_300;
> 
>       --      do 310 k=1,n
>       Do_310 : for K in Weight'range loop -- This hides outer K
>          Weight(K):=Beta(1)*Weight(K)*Weight(K);
>          --  310 continue
>       end loop Do_310;
>       --      return
>       --c
>       --c Set error - no convergence to an eigenvalue after 30
> iterations.
>       --c
>       --  400 ierr=l
>       --      return
>       --      end
>    end Gauss;
> begin -- Make_It_Compile;
>    null;
> end Make_It_Compile;
> 
> --
> Jeff Carter
> "You cheesy lot of second-hand electric donkey-bottom biters."
> Monty Python & the Holy Grail


-- 

Gary Scott
mailto:scottg@flash.net

mailto:webmaster@fortranlib.com
http://www.fortranlib.com

Support the GNU Fortran G95 Project:  http://g95.sourceforge.net

Re: Help with translation: Fortran77 => Ada95

[tmoran]

>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;

Re: Help with translation: Fortran77 => Ada95

[Jeffrey Carter]

Gary Scott wrote:
> 
> Interesting that it's almost exactly twice as long as the Fortran.

Interesting that it compiles and the FORTRAN doesn't. Interesting that,
if you leave off all the stuff needed to make it compile, and all the
comments showing the original FORTRAN, and all the extra white space not
found in the FORTRAN, it's certainly no longer than the FORTRAN.

Try comparing kumquats to quinces next time.

-- 
Jeffrey Carter

Re: Help with translation: Fortran77 => Ada95

[Gary Scott]

Hi,

Jeffrey Carter wrote:
> 
> Gary Scott wrote:
> >
> > Interesting that it's almost exactly twice as long as the Fortran.
> 
> Interesting that it compiles and the FORTRAN doesn't. Interesting that,
> if you leave off all the stuff needed to make it compile, and all the
> comments showing the original FORTRAN, and all the extra white space not
> found in the FORTRAN, it's certainly no longer than the FORTRAN.
> 
> Try comparing kumquats to quinces next time.

Apparently, you took my comments as a personal attack.  I never make
personal attacks, I merely intended to note the difference in style and
a slight tendency to being verbose of the posted update as compared to
the original.  I realized that the original source was poorly designed,
was incomplete, and had several flat-out errors.  The original was not
poorly designed because of the language used, necessarily, and could be
dramatically improved using the same language or a newer variant. 
Perhaps, this was not a point worth making in this particular news
group, but it is just a news group...

> 
> --
> Jeffrey Carter


-- 

Gary Scott
mailto:scottg@flash.net

Re: Help with translation: Fortran77 => Ada95

[Gary Scott]

Hi,

I guess I wasn't intending to start a competition to see how few
statements this could be written in.  I simply noted that the non-blank,
non-comment portion of the posted example was about twice as long as the
original, no big deal.  I also noted that significant improvements could
be made to the original.  Here's a minimalist update which only
eliminates the GOTOs.  At least it's more readable to me, if no one else
(not 100% guaranteed, certainly there's still lots of room for
improvement, a few opportunities for  whole or partial array
operations):

!
!        Input:  n - - the number of points in the Gaussian quadrature
!                      formula; type integer
!                alpha,beta - - arrays of dimension  n  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 of dimension  n  containing the Gaussian
!                      nodes (in increasing order)  zero(k)=x(k), k=1,2,
!                      ...,n
!                weight - array of dimension  n  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, equal to  -1  if  n  is not in
!                      range, and equal to  -2  if one of the beta's is
!                      negative.
!
! The array  e  is needed for working space.
!
subroutine gauss(n, alpha, beta, eps, zero, weight, ierr)

   implicit none

   integer, intent(in)    :: n                    !Automatic extent
   integer, intent(out)   :: ierr                 !Return code
   integer                :: ii, i, j, k, l, m    !Loop counters

   real, intent(in)       :: alpha(n), beta(n), eps
   real, intent(out)      :: zero(n), weight(n)
   real                   :: e(n), b, c, f, g, p, r, s, mml 
!Intermediates and workspaces

   if (n < 1) then  !Invalid extent
      ierr = -1
      return
   end if

   ierr = 0; zero = 0.0; weight = 0.0; e = 0.0 !Initialize

   zero(1) = alpha(1)

   if (beta(1) < 0.0) then
      ierr = -2
      return
   end if

   weight(1) = beta(1)

   if (n == 1) return

   weight(1) = 1.0

   do k = 2, n
      zero(k) = alpha(k)

      if (beta(k) < 0.0) then
        ierr = -2
        return
      end if

      e(k - 1) = sqrt( beta(k) )
      weight(k) = 0.0
   end do

   j = 0

   do l = 1, n

      do m = l, n           ! Look for a small subdiagonal element
         if (m == n) exit
         if (abs( e(m) ) <= eps * (abs( zero(m) ) + abs( zero(m + 1) )))
exit
      end do

      p = zero(l)

      if (m == l) then
         j = 0
         cycle
      end if

      if (j == 30) then !No convergence after 30 interations
         ierr = 1
         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=m-1 step -1 until l do ...
!
      do ii = 1, mml

         i = m - ii
         f = s * e(i)
         b = c * e(i)

         if (abs(f) < abs(g)) then
            s = f / g
            r = sqrt(s * s + 1.0)
            e(i + 1) = g * r
            c = 1.0 / r
            s = s * c
         else
            c = g / f
            r = sqrt(c * c + 1.0)
            e(i + 1) = f * r
            s = 1.0 / r
            c = c * s
         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 do

      zero(l) = zero(l) - p
      e(l) = g
      e(m) = 0.0

   end do
!
! Order eigenvalues and eigenvectors.
!
   do ii = 2, n
      i = ii - 1
      k = i
      p = zero(i)

      do j = ii, n
         if (zero(j) >= p) cycle
         k = j
         p = zero(j)
      end do

      if (k >= i) cycle

      zero(k) = zero(i)
      zero(i) = p
      p = weight(i)
      weight(i) = weight(k)
      weight(k) = p
   end do

   do k = 1, n
      weight(k) = beta(1) * weight(k) * weight(k)
   end do

   return

end subroutine




tmoran@acm.org wrote:
> 
> >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.
> 

<snip>



-- 

Gary Scott
mailto:scottg@flash.net

mailto:webmaster@fortranlib.com
http://www.fortranlib.com

Support the GNU Fortran G95 Project:  http://g95.sourceforge.net