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

Перфокарточный Fortran

Сейчас работаю с программным пакетом, который использует набор библиотек. Одна из них обсчитывает геофизику (волны Лява и Рэлея). Обнаружил, что она сохранилась ещё со времён перфокарт и магнитных лент, будучи написанной в соответствии со стандартными бланками разметки (на исходную публикацию этого пакета ПО в 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

Кажется, пора заводить жанр героического эпоса о борьбе отважных программистов и чудовищ из глубины веков. Программист, приносящий Машине хороший выверенный код, отбитый в схватке с ордами перфокарт, должен занять своё место наравне с рыцарями, отбивающих у драконов прекрасных дам. Герой даже будет чернокнижником и некромантом, что только добавит пикантности сюжету.

Желаю приятного субботнего вечера!

Комментариев нет:

Отправить комментарий