Сейчас работаю с программным пакетом, который использует набор библиотек. Одна из них обсчитывает геофизику (волны Лява и Рэлея). Обнаружил, что она сохранилась ещё со времён перфокарт и магнитных лент, будучи написанной в соответствии со стандартными бланками разметки (на исходную публикацию этого пакета ПО в 1978 году можете посмотреть здесь - в конце более чем двухсотстраничного документа приложена переводная (с русского!) статья по математике, которая призвана помочь разобраться с идеями, заложенными в программе).
И да, с тех пор его так и используют. Я попытался найти более свежую версию, но альтернатив не нашлось. Другие найденные мной программы используют вычислительное ядро в виде этого же оригинального пакета. Кстати, за что я очень ценю фортран: этот древний код даже вперемежку с кодом новых стандартов (а свежайший - это 2008 год) совершенно спокойно компилируют оба основных современных компилятора - gfortran (GCC) и ifort. Покажите мне любой другой язык программирования, для которого поддерживают прямую совместимость фрагментов с сорокалетней разницей в возрасте и стандартах в пределах одного файла. И, кстати, на уровне описания языка поддерживающий тип данных комплексной переменной и встроенные функции, работающие в мнимых числах (что очень актуально для любой вычислительной математики, работающей с волнами).
Тем не менее, поддерживать этот код, а тем более - модицировать (например, поправки на сейсмические затухания), достаточно сложно. Поэтому я решил его переписать на современный стандарт, благо 80 килобайт - это не так много.
Код писан в те времена, когда переходы по строкам были очень дешёвы (для компьютера), а вычисления - очень дороги. В функциях вычисления собственных значений тут до двухста переменных, хранящих косинусы и синусы. А вот циклов и условных переходов в современном понимании нет - есть прыжки по блокам кода. Например:
if (number) 401,402,403
if (number1.lt.number2) goto 505
Если number меньше нуля, код идёт на позицию 401, если равен нулю - то на 402, если же больше - то на 403. Второе выражение выполняет переход на метку 505, если number1 меньше (less than) number2. И ещё:
10 CONTINUE
Означает просто продолжить выполнение (перейдя на следующую строчку). Однако если это метка завершения цикла, то выполнение возвращается на первую строчку:
DO 10 I=1,20
Так задаются циклы для всех i от 1 до 20, причём последняя инструкция цикла находится на строке с меткой 10.
Теперь же я хочу всем желающим проникнуться духом тёплой ламповой эпохи, а так же некрофилам от лингвистики и литературы - перевести на любой известный язык программирования отрывок кода ниже:
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********************************************************************
Вообще, очень интересное (для программирующих) упражнение, которое даёт очень наглядное ощущение того, как далеко ушли компьютеры от машины Тьюринга. Знакомые же с языком ассемблера заметят, насколько к нему близок этот перфокарточный фортран. Отдельного внимания заслуживает навигация по блокам кода для максимального переиспользования фрагментов кода - скорее всего, для минимизации ошибок при копировании. В то время как прыжок в неправильный блок моментально приведёт к заведомо бредовому ответу, пропуск одной инструкции может вести к сложно выявляющимся ошибкам.
Из предыдущего фрагмента у меня получился такой код (он тут уже несколько оптимизирован под конкретную процедуру):
Из предыдущего фрагмента у меня получился такой код (он тут уже несколько оптимизирован под конкретную процедуру):
! 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
Кажется, пора заводить жанр героического эпоса о борьбе отважных программистов и чудовищ из глубины веков. Программист, приносящий Машине хороший выверенный код, отбитый в схватке с ордами перфокарт, должен занять своё место наравне с рыцарями, отбивающих у драконов прекрасных дам. Герой даже будет чернокнижником и некромантом, что только добавит пикантности сюжету.
Желаю приятного субботнего вечера!
Комментариев нет:
Отправить комментарий