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=-0.3 required=5.0 tests=BAYES_00, REPLYTO_WITHOUT_TO_CC autolearn=no 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 18:34:03 PST Path: archiver1.google.com!newsfeed.google.com!newsfeed.stanford.edu!feed.textport.net!news.stealth.net!newscon02.news.prodigy.com!newsmst01.news.prodigy.com!prodigy.com!postmaster.news.prodigy.com!newssvr15.news.prodigy.com.POSTED!not-for-mail Message-ID: <3B43C526.991E0144@flash.net> From: Gary Scott Reply-To: scottg@flash.net Organization: Home X-Mailer: Mozilla 4.77 [en]C-DIAL (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> <3B435888.55E1BE36@acm.org> Content-Type: text/plain; charset=us-ascii Content-Transfer-Encoding: 7bit NNTP-Posting-Host: 64.48.221.221 X-Complaints-To: abuse@prodigy.net X-Trace: newssvr15.news.prodigy.com 994296782 6242077 64.48.221.221 (Wed, 04 Jul 2001 21:33:02 EDT) NNTP-Posting-Date: Wed, 04 Jul 2001 21:33:02 EDT Date: Thu, 05 Jul 2001 01:33:02 GMT Xref: archiver1.google.com comp.lang.ada:9447 Date: 2001-07-05T01:33:02+00:00 List-Id: 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