From mboxrd@z Thu Jan 1 00:00:00 1970 X-Spam-Checker-Version: SpamAssassin 3.4.4 (2020-01-24) on polar.synack.me X-Spam-Level: X-Spam-Status: No, score=-1.9 required=5.0 tests=BAYES_00 autolearn=ham autolearn_force=no version=3.4.4 X-Google-Language: ENGLISH,ASCII-7-bit X-Google-Thread: 103376,83d867c7a987cc31 X-Google-Attributes: gid103376,public X-Google-ArrivalTime: 2001-07-04 10:55:30 PST Path: archiver1.google.com!newsfeed.google.com!newsfeed.stanford.edu!nntp.cs.ubc.ca!newsfeed.direct.ca!look.ca!newsfeed1.earthlink.net!newsfeed.earthlink.net!newsmaster1.prod.itd.earthlink.net!newsread2.prod.itd.earthlink.net.POSTED!not-for-mail Message-ID: <3B435888.55E1BE36@acm.org> From: Jeffrey Carter X-Mailer: Mozilla 4.7 [en] (Win98; U) X-Accept-Language: en MIME-Version: 1.0 Newsgroups: comp.lang.ada Subject: Re: Help with translation: Fortran77 => Ada95 References: <9hun38$t0m$1@ulysses.noc.ntua.gr> Content-Type: text/plain; charset=us-ascii Content-Transfer-Encoding: 7bit Date: Wed, 04 Jul 2001 17:55:34 GMT NNTP-Posting-Host: 206.133.139.107 X-Complaints-To: abuse@earthlink.net X-Trace: newsread2.prod.itd.earthlink.net 994269334 206.133.139.107 (Wed, 04 Jul 2001 10:55:34 PDT) NNTP-Posting-Date: Wed, 04 Jul 2001 10:55:34 PDT Organization: EarthLink Inc. -- http://www.EarthLink.net X-Received-Date: Wed, 04 Jul 2001 10:53:12 PDT (newsmaster1.prod.itd.earthlink.net) Xref: archiver1.google.com comp.lang.ada:9440 Date: 2001-07-04T17:55:34+00:00 List-Id: 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