## суббота, 21 апреля 2018 г.

### Fortran from punch card era

I was working on one geophysical software, when I found one of the libraries to compute Love and Rayleigh waves is still from the times of punch cards and tapes. Its source fits the standard of fortran pads (here you may take a look on the original release of this utility in the year 1978, including "a translation of a Russian article" to "help understand the theory used").

And yes, this piece of code is still used. I tried to look for another library, but haven't found any (at least, in public domain). All the other source codes I've seen use the same tool. That's actually why do I appreciate Fortran so much: please, show me any other programming language, so that both widely used compilers (GCC and ifort) can assemble a working program from a file consisting a mixture of standards and codes written with 40 years difference! And also having a native datatype of complex variable and functions of is, which is very needed by any wave theory math.

However, it's difficult to manage this code and improve it (for example, to add seismic wave attenuation corrections etc). Therefore, I decided to rewrite it to modern day standard. 80 kilobytes is not so lot.

This code was written when jumps between code pieces was very cheap relatively to the actual computations of some value. Eigenfunctions in this code have up to 200 variables to store sines and cosines. However, there no conditional cases or cycles in our modern-day meaning. There are jumps to certain labels of positions in code. For example:
if (number) 401,402,403
if (number1.lt.number2) goto 505
If number is less than zero, code goes to label 401, if it is zero, it goes to 402, and to 403 in the rest of values. The second expression jumps to label 505 is number1 is less than number2. And also:
10       CONTINUE
Means to continue execution (go to the next line). In the case of cycle it return to the first line of it:
DO 10 I=1,20
Do it for all the i-s in between 1 and 20, and the last statement is on the line with label 10.
Now I'd like to all the people, as well as linguistics/literature necrophiliacs, who want to feel the spirit of the good old days to translate this code to any modern language:
C********************************************************************
C---------------start of loop for Period k--------------------------S
do 9998 k=1,kmax
call flat1(d,rho,a,b,mmax,kind)
C--------------------loop for modes---------------------------------S
do 9997 iq=1,kmode
if(k-1) 605,605,599
599     if(iq-2) 600,601,601
600     c1=0.90*c(k-1,1)
goto 605
601     if(dc*(c(k-1,iq)-c(k,iq-1))) 602,603,604
602     c1=c(k,iq-1)+0.01*dc
goto 605
603     c1=c(k,iq-1)+0.01*dc
goto 605
604     c1=c(k-1,iq)
605     continue
idrop=0
999     del1=dltar(c1,t1,ifunc)
80      c2=c1+dc
idrop=0
del2=dltar(c2,t1,ifunc)
if(sign(1.,del1).ne.sign(1.,del2)) goto 54
81      c1=c2
del1=del2
qbq=0.8*b(1)
if(c1-0.8*b(1)) 250,251,251
251     if(c1-(b(mmax)+0.3)) 252,250,250
252     goto 80
54      continue
call nevill (t1,c1,c2,del1,del2,ifunc,cn)
if(lstop.eq.0)                  go to 3006
lstop=0
if(k.eq.1.and.iq.eq.1)          go to 3005
if(k.eq.1)                      go to 3004
do 3003 i=1,k-1
do 3000 j=1,mode
if(i.gt.imax(j))                go to 3001
3000    continue
jmax=mode
go to 3002
3001    jmax=j-1
3002    continue
3003    continue
3004    if(iq.eq.1)                     go to 3005
jmax=iq-1
3005    continue
go to 9999    ! to the end of the subroutine
3006    c1=cn
if(c1-b(mmax)) 121,121,250
121     c(k,iq)=c1
C----------------ifunc=kind : Rayleigh if =2 or Love if =1-----------
goto (6001,6002),ifunc
6002    continue
ratio(k,iq)=dltar(c1,t1,3)
goto 6003
6001    continue
6003    continue
c1=c1+0.01*dc
imax(iq)=k
goto 9997
250     if(k*iq-1) 256,256,255
256     continue
write(*,*) ncount
if(ncount.eq.1)then
PRINT 258
do i=1,mmax
print *, b(i), a(i) ,rho(i), d(i), qs_ref(i)
enddo
258     format(1x,/22x,'Improper initial value. No zero found.')
endif
goto 9999  ! this goes to the end of the subroutine
9997    continue
C-----------------end of loop for mode------------------------------E
C--------------------search of root is finished----------------------
goto 9998
255     kmode=iq-1
if(kmode) 9996,9996,9998
9998    continue
C-----------------end of loop for period k--------------------------E
9996    continue
C********************************************************************

That's a pretty interesting exercise which gives a very good feeling of the long way made by computers from Turing's machine to the modern day gadgets. Persons who have coded in assembly language may also see how close is this code to it. You may also notice the concept of re-use: they tried to reuse any single piece of code. That makes sense: while one can make typos in one of the lines, leading to a hard-to-find bugs (like -0.001 instead of +0.001), usage of a wrong block will be clearly seen because of weird results.
For the aforementioned piece of code I go this listing (a bit optimised for the routine is uses):
! start of loop for Period k
do k=1,kmax
! kmode is = mode = 1 in our case (fundamental modes)
do iq=1,kmode
if (k-1 > 0) then
if (iq-2 < 0) then
c1=0.90*cc(k-1,1)
else
if(dc*(cc(k-1,iq)-cc(k,iq-1)) <= 0.) then
c1=cc(k,iq-1)+0.01*dc
else
c1=cc(k-1,iq)
endif
endif
endif
idrop=0 ; del1=dltar(c1,t1,ifunc)
goin = .false.
do
c2=c1+dc
idrop=0 ; del2=dltar(c2,t1,ifunc)
if (sign(1., real(del1))*sign(1., real(del2)) < 0.) then
call nevill (t1,c1,c2,del1,del2,ifunc,cn)
goin = .true.
exit
endif
c1=c2
del1=del2
qbq=0.8*b(1)
if (c1-0.8*b(1) < 0. .and. c1-(b(mmax)+0.3) >= 0.) exit
enddo
if (goin) then
if (lstop /= 0) then
lstop=0
if(k /= 1) then
do i=1,k-1
jmax=mode
do j=1,mode
if(i > imax(j)) then
jmax=j-1
exit
endif
enddo
enddo
endif
if(iq.eq.1) return
jmax=iq-1
endif
c1 = cn
if(c1-b(mmax) <= 0) then
cc(k,iq)=c1
if (ifunc == 2) ratio(k,iq)=dltar(c1,t1,3) ! 2 = Ray, 1 = Love
c1=c1+0.01*dc
imax(iq)=k
cycle
endif
endif
if(k*iq-1>0) exit
write(*,*) ncount
if(ncount == 1)then
PRINT *, 'Improper initial value. No zero found in calcul.f'
do i=1,mmax
print *, b(i), a(i) ,rho(i), d(i), qs_ref(i)
enddo
endif
return
enddo
enddo ! end of loop for k


I think humankind needs a new wave of epic poetry. A programmer will bring to the Machine a list of code captured in a fight with with hordes of wild punch cards - similar to a knight bringing a head of defeated dragon. The hero will be also a necromancer and practice sorcery.

And yeah, have a nice Saturday evening!

#### 1 комментарий:

1. 