From: Gary Scott <scottg@flash.net>
Subject: Re: Help with translation: Fortran77 => Ada95
Date: Thu, 05 Jul 2001 01:33:02 GMT
Date: 2001-07-05T01:33:02+00:00 [thread overview]
Message-ID: <3B43C526.991E0144@flash.net> (raw)
In-Reply-To: 3B435888.55E1BE36@acm.org
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
next prev parent reply other threads:[~2001-07-05 1:33 UTC|newest]
Thread overview: 11+ messages / expand[flat|nested] mbox.gz Atom feed top
2001-07-04 9:19 Help with translation: Fortran77 => Ada95 N&J
2001-07-04 17:04 ` Warren W. Gay VE3WWG
2001-07-04 17:17 ` Colin Paul Gloster
2001-07-04 17:55 ` Jeffrey Carter
2001-07-05 1:33 ` Gary Scott [this message]
2001-07-05 6:03 ` tmoran
2001-07-06 1:35 ` Gary Scott
2001-07-05 17:36 ` Jeffrey Carter
2001-07-06 1:14 ` Gary Scott
2001-07-04 19:51 ` [SOLVED] " N&J
2001-07-04 23:25 ` GianLuigi Piacentini
replies disabled
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox