Fun with Fortran on the PiDP-11

[Peter campbell-burns]

I last wrote programs in Fortran in the early 1980s - the applicatiion was digitisation of bubble chamber photographs and the platform was PDP-11 / RSX-11M/Plus.  To squeeze a quart into the proverbial pint pot there was heavy use of overlays - something that today's coders don't need to be concerned about.  I remember helpful messages from the task builder "Error in file ." - whatever "." was.  My employer was University College London; it was without question the best job I ever had, but also the most poorly paid given the expense of living in our capital city.

Some time ago I had a go at modelling dark flight of meteorites; my interest here is that I am a member and co-founder of the UK Meteor Network (UKMON).  We hoped one day to recover a fresh meteorite fall.  This happened last year with the Winchcombe meteorite - but alas not using my amateurish dark flight code.  My code was written using object-oriented Delphi/Pascal.  Porting this to my PiDP-11 sounded like a fun project.

I approached the project in two stages - first I ported the code to gfortran (only because I had a linux laptop whichmeant I could work on the move) and then to DEC Fortran77 on my PiDP-11.   Expected a lot of differences between DEC and gfortran flavours of fortran 77 and was not disappointed.  It was an interesting experience.   Code that compiles on the gfortran flavour of F77 required a lot of further changes. 

This exercise has done a lot to remind reminded me just how far things have changed and the constraints we used to work within - not just limitations on memory and storage but even seemingly simple things like length variable names.  My Delphi solution read metological data into a large array; not something I want to do on a PDP-11.

Well, it compiles, it runs, but is so far producing incorrect results - I probably made a mistake somewhere shortening variable names of changing logical statements.  The Fortran fun continues - this will keep me busy for a little while longer. 

[Johnny Billquist]

Two things. Remember to use split I/D space, which makes a big
difference for when you start needing overlays. And use supervisor mode
libraries, which also removes a lot of your memory pressure.

Also, if you are dealing with large arrays in FORTRAN, you might want to
look at virtual arrays... ;-)

Apart from that, sounds like you are having fun. Yes, the 6 character
limitation on global symbols (and local for FORTRAN) can be a bit limiting.

Anyway, if you have questions, just ask.

Johnny

> --
> You received this message because you are subscribed to the Google
> Groups "[PiDP-11]" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to pidp-11+u...@googlegroups.com
> <mailto:pidp-11+u...@googlegroups.com>.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/pidp-11/458c38a5-584e-4280-b239-87ec9fc2d933n%40googlegroups.com
> <https://groups.google.com/d/msgid/pidp-11/458c38a5-584e-4280-b239-87ec9fc2d933n%40googlegroups.com?utm_medium=email&utm_source=footer>.

--
Johnny Billquist || "I'm on a bus
|| on a psychedelic trip
email: b...@softjar.se || Reading murder books
pdp is alive! || tryin' to stay hip" - B. Idol

[Clem Cole]

Some thoughts below that might make your project go a little more smoothly... [note I do work for Intel in my current day job, and I did work for DEC years before, but these comments are my own not necessarily those of my employer then or now].

On Fri, Mar 4, 2022 at 3:50 AM Peter campbell-burns <pcampbe...@gmail.com> wrote:

I last wrote programs in Fortran in the early 1980s - the applicatiion was digitisation of bubble chamber photographs and the platform was PDP-11 / RSX-11M/Plus.  To squeeze a quart into the proverbial pint pot there was heavy use of overlays - something that today's coders don't need to be concerned about.  I remember helpful messages from the task builder "Error in file ." - whatever "." was.  My employer was University College London; it was without question the best job I ever had, but also the most poorly paid given the expense of living in our capital city.

Some time ago I had a go at modelling dark flight of meteorites; my interest here is that I am a member and co-founder of the UK Meteor Network (UKMON).  We hoped one day to recover a fresh meteorite fall.  This happened last year with the Winchcombe meteorite - but alas not using my amateurish dark flight code.  My code was written using object-oriented Delphi/Pascal.  Porting this to my PiDP-11 sounded like a fun project.

Do you know about the FreePascal distribution, which accepts Delphi as one of its many input languages?   While FreePascal runs on Linux, Mac, Windows, the old Omsi Pascal compiler runs on the PDP-11 and associated OS's and binaries of the same are still floating around the Internet these days  (I leave it to you to do your own search).  But, you might try tweaking the Delphi code to be closer to Pascal code to make it more acceptable to the Omsi compiler as one solution, using the FreePascal compiler to help with that transition.

I approached the project in two stages - first I ported the code to gfortran (only because I had a linux laptop whichmeant I could work on the move) and then to DEC Fortran77 on my PiDP-11.   Expected a lot of differences between DEC and gfortran flavours of fortran 77 and was not disappointed.  It was an interesting experience.   Code that compiles on the gfortran flavour of F77 required a lot of further changes. 

Ah, a possible tactical error here, trying to use gfortran as the intermediate step could be loaded with many difficulties/hidden dragons.

It turns out that the old DEC compiler front-end and associated runtimes still live and are maintained/continue to be enhanced.  Remember Intel ended up with the DEC compiler assets and people and their compiler have been using the exact same front-end as was originally written years ago [little known factoid from my friends in the compiler team, the DEC FTN Frontend was/is written in Pascal].  The target IL and code generators have changed, but the front-end and runtime are still accepting 'dust decks' just like they did years ago on the DEC systems.

That said, the current back end for it is LLVM based and the tool is now called ifx (ifort, the previous Intel Fortran compiler -- a step child of DEC GEM -- can still be found and is still the prefered 'production' compiler at most HPC sites).  Developers can download the new code base for free for on Windows, Linux, Mac (premium support costs $s). It is downloaded as part of the Intel OneAPI HPC Toolkit [you'll also get the traditional debuggers and tuning tools, C, C++, DPC++, HPC Python etc..)

What I would suggest is getting the code running with the new Intel compiler and then slowly making it more DEC FTN-77 like while running on a large address machine with better debug tools.  The current Intel compiler will even accept fixed format FORTRAN-IV [1] as well as modern Fortran-2018 syntax. As you back it up into F77 or F-IV syntax using the Intel compiler, later trying it with the original DEC Compiler the differences should start to be more ISA based issues and I would expect them to be easier to identify.

1] See my favorite test - which we call the 'Eklund test.'  Dave was the fixed format guru in DEC's compiler team until he retired from Intel a couple of years ago.  This was posted to another mailing list a few years back, so I submitted it in Dave honor to the FTN team:

C    This FORTRAN program may be compiled and run on a Norsk Data
C    computer running SINTRAN and the FTN compiler.  It uses only
C    FORTRAN reserved words, and contains just one numerical
C    constant, in a character string (a format specifier).  When
C    you run it, it prints a well known mathematical construct...
C
C    Even FORTRAN is a block structured programming language:
C
      PROGRAM
     ;PROGRAM;INTEGERIF,INTEGER,GOTO,IMPLICIT;REALREAL,DIMENSION,EXTERNA
     AL,FORMAT,END;INTEGERLOGICAL;REALCOMPLEX,DATA,CALL,ASSIGN,CHARACTER
     R;DOFORIF=INTEGER,INTEGER;ENDDO;INTEGER=IF+IF;GOTO=INTEGER*INTEGER*
     *INTEGER*INTEGER-INTEGER-IF;CALLFUNCTION(IMPLICIT,REAL,DIMENSION,EX
     XTERNAL,FORMAT,END,LOGICAL,COMPLEX,DATA,CALL,ASSIGN,CHARACTER);CALL
     LSUBROUTINE(IMPLICIT,LOGICAL,GOTO,IF,INTEGER)
      END
      SUBROUTINEFUNCTIO
     ON(IMPLICIT,REAL,DIMENSION,EXTERNAL,FORMAT,END,LOGICAL,COMPLEX,DATA
     A,CALL,ASSIGN,CHARACTER);RETURN
      END
      SUBROUTINESUBROUTINE(IMPLICIT,L
     LOGICAL,GOTO,IF,INTEGER);INTEGERGOTO,IMPLICIT(GOTO),LOGICAL(GOTO),I
     IF,INTEGER,EXTERNAL,RETURN;DOFOREXTERNAL=IF,GOTO;DOFORRETURN=INTEGE
     ER,EXTERNAL-IF;IMPLICIT(RETURN)=LOGICAL(RETURN)+LOGICAL(RETURN-IF);
     ;ENDDO;IMPLICIT(IF)=IF;IMPLICIT(EXTERNAL)=IF;DOFORRETURN=IF,GOTO-EX
     XTERNAL;WRITE(IF,'(''$  '')');ENDDO;DOFORRETURN=IF,EXTERNAL;WRITE(I
     IF,'(''$''I4)')IMPLICIT(RETURN);ENDDO;WRITE(IF,'( /)');DOFORRETURN=
     =IF,GOTO;LOGICAL(RETURN)=IMPLICIT(RETURN);ENDDO;ENDDO
      END

[johntk...@gmail.com]

Writing Fortran apps sounds like fun - and a lot less fighting the OS that modern devs have to do (although overlays sound a headache). 

BTW if tracking meteors is a curve fitting problem, there are some great machine learning models that would be a good tool to try. 

[Johnny Billquist]

FORTRAN isn't that bad, and DEC's dialects are among the nicer, if you
ask me.

Overlays can be a headache, but not necessarily. It's all about
describing the dependencies between modules. The larger your software
is, the more complex that can become. But in many cases, it can be very
straight forward.

Johnny

> --
> You received this message because you are subscribed to the Google
> Groups "[PiDP-11]" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to pidp-11+u...@googlegroups.com
> <mailto:pidp-11+u...@googlegroups.com>.
> To view this discussion on the web visit

> https://groups.google.com/d/msgid/pidp-11/aca8cfad-f962-4170-a1bd-115e8b35e9d7n%40googlegroups.com
> <https://groups.google.com/d/msgid/pidp-11/aca8cfad-f962-4170-a1bd-115e8b35e9d7n%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Neal G.]

About Virtual Arrays ... is there anything beyond replacing DIMENSION with VIRTUAL that needs to be done to use them?

The docs seem to indicate that this is all that is necessary; yet I encountered errors in a small test application.

The application runs fine when using a normal array, but encounters this error when using a virtual array.

TT3  --  ERROR 112
Virtual array mapping error
  in   "HSORT " at or after 24
  from "UHEAP4" at or after 13

TT3  --  Exiting due to ERROR 3
Odd address trap (SST0)
at PC = 003070
  in   "HSORT " at or after 24
  from "UHEAP4" at or after 13

- Neal G.

[Johnny Billquist]

I think you need to open a file that is the backing store, and associate it.

But I don't remember all the details offhand. Why not just check the manual?

Johnny

> <https://groups.google.com/d/msgid/pidp-11/aca8cfad-f962-4170-a1bd-115e8b35e9d7n%40googlegroups.com?utm_medium=email&utm_source=footer

> <https://groups.google.com/d/msgid/pidp-11/aca8cfad-f962-4170-a1bd-115e8b35e9d7n%40googlegroups.com?utm_medium=email&utm_source=footer>>.
>
>
> --
> Johnny Billquist || "I'm on a bus
> || on a psychedelic trip
> email: b...@softjar.se || Reading murder books
> pdp is alive! || tryin' to stay hip" - B. Idol
>

> --
> You received this message because you are subscribed to the Google
> Groups "[PiDP-11]" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to pidp-11+u...@googlegroups.com
> <mailto:pidp-11+u...@googlegroups.com>.
> To view this discussion on the web visit

> https://groups.google.com/d/msgid/pidp-11/59d64292-27a9-40d1-a4b5-47851cf47557n%40googlegroups.com
> <https://groups.google.com/d/msgid/pidp-11/59d64292-27a9-40d1-a4b5-47851cf47557n%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Johnny Billquist]

Sorry. I was thinking of virtual arrays in BASIC+2. Seems they would be
handled without a file backing in F77.

I don't know why you are getting that error. Is the task build
checkpointable?

Johnny

[Clem Cole]

On Wed, Mar 9, 2022 at 12:19 PM Neal G. <cven...@earthlink.net> wrote:

About Virtual Arrays ... is there anything beyond replacing DIMENSION with VIRTUAL that needs to be done to use them?

The docs seem to indicate that this is all that is necessary; yet I encountered errors in a small test application.

The application runs fine when using a normal array, but encounters this error when using a virtual array.

TT3  --  ERROR 112
Virtual array mapping error
  in   "HSORT " at or after 24
  from "UHEAP4" at or after 13

TT3  --  Exiting due to ERROR 3
Odd address trap (SST0)
at PC = 003070
  in   "HSORT " at or after 24
  from "UHEAP4" at or after 13

- Neal G.

Which version of the compiler are you using?

http://www.bitsavers.org/pdf/dec/pdp11/rt11/v4.0_Mar80/4_fortran/DEC-11-LFLRA_FORTRAN_Language_Reference_Manual_Jun77.pdf

Appendix D describes Virtual arrays.    It's been years (late 1970s) since I played with these and then it was on RT-11 (not RSX) - although the compilers were supposed to be mostly the same (I was told by folks in TLG that they were built from the same BLISS-11 sources at one point).  And yes, to declare them you replace the DIMENSION statement with VIRTUAL --- but there are rules which I have completely forgotten WRT to size, and you can't use them with block COMMON and other such -  Check the manual and good luck. You might be violating one of the restrictions.   As I said, I have completely forgotten what we did to make it all work (we were developing some early graphics dithering algorithms and needed some large arrays). The code had started on a CDC box and moved back to the PDP-11.   I do remember the code finally got moved to a Vax later.

[Mark Matlock]

Neal,

   I just checked an old program I wrote that used F77 virtual arrays. Just to make sure nothing had changed much I recompiled it and rebuilt it and ran it and everything ran fine. The errors you are seeing might result from out of bounds array references. I noticed in my old build command file that it accepted a command line debug switch (e.g. @CALIB.CMD D ) which compiled the F77 with array bounds checking “/CK” that generates extra code to test that each attempt to reference an array element is valid. It also links in the PDP-11 Symbolic debugger and turns off compiler optimization.

   I had some additional programs and data files (this was all part of my Master’s thesis) that ran on a PDP-11/44 in the day. It lacked floating point and we ran a DECUS FPP emulator so our floating point was slow. A batch file that ran all the data through the suite of programs took a little over an hour. On my 11/83 today it takes 26 minutes and on a Mentec M100-04 (~11/93) it takes 19 minutes. Their were a number of generated FTIR spectrum plots that were spooled to a HP7550 plotter from the batch job.

Best,

Mark

Here is the F77 code:

       PROGRAM CALIB
      PARAMETER (MS=15, MP=1000, MC=4)
      VIRTUAL R(MC+1,MP), S(MS,MP), U(MP), W(MP)
      VIRTUAL X(MC+1,MS), Y(MS), Z(MC+1,MS)
      REAL*8 R, S, U, X, Y, Z
      REAL*8 DET, SE, YD, YE, C(MS,MC), T(MC+1)
      INTEGER*2 I, J, K, NS, NP, NC
      CHARACTER FNAM*24, TITLE*40
C
C     NS Number of Spectra
C     MS Maximum Number of Spectra
C     NP Number of Points in a Spectrum
C     MP Maximum Number of Points in a Spectrum
C     NC Number of Components in the Sample Composition
C     MC Maximum Number of Components
C     S  Spectra Matrix     NS spectra x NP points
C     C  Composition Matrix NS spectra x NC components
C     X  Regression Matrix  NS spectra x NC+1 components
C     Y  Regression Matrix  NS spectra x a spectral point
C     R  Regression Coefficient Matrix NC+1 x NP
C     T  Temporary Coefficient Storage
C     U  Uncertainty Vector NP values = Standard Error
C        of Estimate for each Regression performed
C
      TYPE 100,'Enter Name of parameter file '
  100 FORMAT (1X,A,$)
      ACCEPT 110,FNAM
  110 FORMAT(A)
      OPEN (UNIT=2, STATUS='OLD', FILE=FNAM,
     +CARRIAGECONTROL='LIST', READONLY)
C
C     The parameter file has the Number of Spectra,
C     the Number of points in each Spectrum, and
C     the number of components, the each record
C     after that is the name of one of the NS files
C     containing the Spectral Data
C
      READ (2,*) NS, NP, NC
      DO 140 I=1,NS
      READ (2,120) FNAM
  120 FORMAT (A)
C
      TYPE *,'Opening File ',FNAM
      OPEN (UNIT=1, STATUS='OLD', FILE=FNAM,
     +CARRIAGECONTROL='LIST', READONLY)
C
C     For each Spectrum read the Title,
C     the Component Levels and the NP point spectrum
C
      READ (1,120) TITLE
      READ (1,*) (C(I,J),J=1,NC)
      DO 130 J=1,NP
  130 READ (1,*) W(J), S(I,J)
      TYPE *,'Finished reading Data for ',TITLE
  140 CLOSE (UNIT=1)
      CLOSE (UNIT=2)
C
C     Build the Component Matrix X used for all Regressions
C
      DO 150 I=1,NS
      X(1,I)=1.
      DO 150 J=1,NC
  150 X(J+1,I)=C(I,J)
C
C     Build the Y Matrix for each spectral point
C     And perform a regression
C
      DO 190 I=1,NP
      DO 160 J=1,NS
  160 Y(J)=S(J,I)
C
      CALL REGRES (X, Y, NS, NC+1, T, DET, Z)
      IF (DET.EQ.0.0) THEN
      TYPE 170
  170 FORMAT (' Regression Failed',/,
     +' Sample Compositions are Colinear')
      STOP
      ENDIF
      DO 180 J=1,NC+1
  180 R(J,I)=T(J)
  190 CONTINUE
      TYPE *,'Determinant is ',DET
C
C     Calculate the Uncertainties which are the RMS averages
C     of the deviations of the predictions from the actual
C     spectra, also known as the standard error of estimate
C
      DO 220 I=1,NP
      SE=0.
      DO 210 J=1,NS
      YE=0.
      DO 200 K=1,NC+1
  200 YE=YE+X(K,J)*R(K,I)
      YD=S(J,I)-YE
  210 SE=SE+YD*YD
  220 U(I)=SQRT(ABS(SE/NS))
C
C     Store the Regression Coefficients in a file that has
C     NC+1 coefficients for each record plus an uncertainty
C     and NS records for each point in the spectra
C
      OPEN (UNIT=1, STATUS='NEW', FILE='REGRES.DAT',
     +CARRIAGECONTROL='LIST')
C
      WRITE (1,230) NP, NC
  230 FORMAT (2I6)
  240 FORMAT (F8.0,5F12.6)
      DO 250 I=1,NP
  250 WRITE (1,240) W(I),(R(J,I),J=1,NC+1),U(I)
      CLOSE (UNIT=1)
      END

Here is the build command file:

.IF P1 EQ "D" F77 CALIB/DS=CALIB/-OP/DB/CK

.IF P1 NE "D" F77 CALIB/DS=CALIB

.IF <EXSTAT> <> 1 .EXIT

.IF P1 EQ "D" F77 REGRES/DS=REGRES/-OP/DB/CK

.IF P1 NE "D" F77 REGRES/DS=REGRES

.IF <EXSTAT> <> 1 .EXIT

.IF P1 EQ "D" F77 MATINV/DS=MATINV/-OP/DB/CK

.IF P1 NE "D" F77 MATINV/DS=MATINV

.IF <EXSTAT> <> 1 .EXIT

.OPEN CALIB.TKB

.IF P1 EQ "D" .DATA CALIB/ID,,CALIB=CALIB,REGRES,MATINV,LB:[1,1]F77DBG/DA,LB:[1,1]F77FCS/LB

.IF P1 NE "D" .DATA CALIB/ID=CALIB,REGRES,MATINV,LB:[1,1]F77FCS/LB

.DATA /

.DATA SUPLIB=FCSFSL:SV

.DATA //

.CLOSE

TKB @CALIB.TKB

PUR *.OBJ,*.TSK,*.STB,*.TKB

[Neal G.]

Thanks everyone for your ideas, information, and examples.

>> Which version of the compiler are you using?

$ FORTRAN/ID
F77 -- PDP-11 FORTRAN-77 V5.4-26

>> I noticed in my old build command file that it accepted a command line

>> debug switch (e.g. @CALIB.CMD D ) which compiled the F77 with array

>> bounds checking “/CK”

The DCR equivalent is "/CHECK". Good suggestion, when compiled with that option I see:

TT3  --  ERROR 94
Array reference outside array

  in   "HSORT " at or after 24
  from "UHEAP4" at or after 13

TT3  --  ERROR 112
Virtual array mapping error
  in   "HSORT " at or after 24
  from "UHEAP4" at or after 13

...

The code I am working with is an adaptation of a heap sort example using an

algorithm which I think originally came from Numerical Recipes for F77. I had

recast it to be compile-able with FORTRAN-IV. Then brought the example to

RSX and RT11 to tinker with F77 in these environments. It looks like something

has become broken in the process. A good time to learn about using the debugger.

>> FTIR spectrum plots

That brings back memories, but not so long ago for me. OMNIC software on Windows/XP,

using a Nicolet 6700 bench. I had a small role on the Engineering team for the 6700 project.

- Neal G.

[Mark Matlock]

Neal,

   If you have downloaded the DU1: DECUS disk from RSX11M.com, look in [6,30]HPSORT.FTN and [6,31]XHPSORT.FTN

They compile ok except for the XHPSORT uses the same 7 character name in the PROGRAM statement and has to be shortened.

   Back in 1983, I was also using a Nicolet FTIR but it had a 20 bit word length Nicolet computer. I used a serial line to transfer spectra

to the PDP-11/44 where it was easier to write my own routines.

Best,

Mark

[Mark Matlock]

For anyone who is interested in the heap sort algorithm, here is some F77 code that demonstrates it, with a timer. It is from

Numerical Recipes for F77 (a great book!). Timing is done with the F77 SECNDS routine which is only accurate to the 60 (or 50) Hz tick.

      PROGRAM xhpsor
C     driver for routine hpsort
      INTEGER i,j
      REAL a(100), t
      open(1,file='TARRAY.DAT',status='OLD')
      read(1,*) (a(i),i=1,100)
      close(1)
C     print original array
      write(*,*) 'Original array:'
      do 11 i=1,10
        write(*,'(1x,10f7.2)') (a(10*(i-1)+j),j=1,10)
11    continue
C     sort array
      t=secnds(0.)
      call hpsort(100,a)
      t=secnds(t)
C     print sorted array
      write(*,*) 'Sorted array:'
      do 12 i=1,10
        write(*,'(1x,10f7.2)') (a(10*(i-1)+j),j=1,10)
12    continue
      write(*,100) t
100   format (' Elapsed Sort Time ',F5.2,' Seconds')
      end

Here is the heap sort routine itself, you can test other sort algorithms by calling your sort routine instead:

      SUBROUTINE hpsort(n,ra)
      INTEGER n
      REAL ra(n)
      INTEGER i,ir,j,l
      REAL rra
      if (n.lt.2) return
      l=n/2+1
      ir=n
10    continue
        if(l.gt.1)then
          l=l-1
          rra=ra(l)
        else
          rra=ra(ir)
          ra(ir)=ra(1)
          ir=ir-1
          if(ir.eq.1)then
            ra(1)=rra
            return
          endif
        endif
        i=l
        j=l+l
20      if(j.le.ir)then
          if(j.lt.ir)then
            if(ra(j).lt.ra(j+1))j=j+1
          endif
          if(rra.lt.ra(j))then
            ra(i)=ra(j)
            i=j
            j=j+j
          else
            j=ir+1
          endif
        goto 20
        endif
        ra(i)=rra
      goto 10

      end

To create random test data set files for sort routine testing. It uses the F77 RAN function which needs an Interger seed

which I get from the SECNDS routine (seconds since midnight):

      PROGRAM tarray
C     create TARRAY.DAT for routine sort routines
      INTEGER i,j
      INTEGER *4 seed
      REAL a(100)
      seed=INT(secnds(0.))
      do 10 i=1,100
10     a(i)=ran(seed)
      open(1,file='TARRAY.DAT',status='NEW')
      write(1,*) (a(i),i=1,100)
      close(1)
C     print array
      write(*,*) 'Random array:'
      do 11 i=1,10
        write(*,'(1x,10f7.2)') (a(10*(i-1)+j),j=1,10)
11    continue
      end

[Neal G.]

Mark,

Yes, [6,30]HPSORT.FTN is the algorithm I had used, revised to sort an integer array, and to replace the IF blocks with GOTO so I can compile with older versions of FORTRAN.

Similar to your example, my test program measures the sort time as a crude performance gauge between different systems, and a comparison with an alternate sort algorithm I had used many years ago.

I swapped the HPSORT.FTN code into my test program and confirmed the problem is present in their code and unrelated to my revisions.

The problem is an unstated limitation, a deficiency, in that Numerical Recipes algorithm. The maximum size of the array the algorithm can sort is 1/2 the maximum value of INTEGER. If this limit is exceeded, one of the index variables will overflow and appear to be negative on the next pass through the algorithm. The algorithm does not detect this out of bounds value and attempts to use the negative index to access the array.

In my test program the array size was 20000. Larger than 1/2 the max value of INTEGER*2. The default DEC F77 INTEGER. So when run, this problem occurs and the memory mapping mechanism which implements the DEC F77 Virtual Array caught this invalid memory access and terminated the program.

When compiled with /CK (DCL /CHECK) the out of bounds situation is detected and reported clearly.

Oddly, when I had used the test program with normal arrays, no error was detected; even which compiled with /CK.

Curious, if the sort was actually working I put the sorted values into a spreadsheet and plotted the result. Yes, even though the sort appeared to complete without error, about one quarter of the values were not sorted correctly.

I'm not entirely certain of the moral of the story; be wary of example code, maybe. But it is certainly enlightening to dive into some of the basic issues developers faced on these classic systems.

- Neal G.

[Rob Pike]

Indeed. Be careful about the code Numerical Recipes if you are doing
important numerical work. My numerical colleagues warned me off it
years ago, as the algorithms there do not worry enough about accuracy
and stability. It would not be described by them as a "great book".

For simple stuff it's fine, of course, but floating point is a
minefield that requires expert navigation.

-rob

> --
> You received this message because you are subscribed to the Google Groups "[PiDP-11]" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to pidp-11+u...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/pidp-11/6eda2dc2-9138-4dbc-8490-09c5c2f648f4n%40googlegroups.com.

[Johnny Billquist]

I'm really surprised that /CK would not have caught this for a normal
array but for a virtual array.

Are you sure about that? Could you provide all the sources? I'm curious
about this.

Johnny

> --
> You received this message because you are subscribed to the Google
> Groups "[PiDP-11]" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to pidp-11+u...@googlegroups.com
> <mailto:pidp-11+u...@googlegroups.com>.
> To view this discussion on the web visit

> https://groups.google.com/d/msgid/pidp-11/6eda2dc2-9138-4dbc-8490-09c5c2f648f4n%40googlegroups.com
> <https://groups.google.com/d/msgid/pidp-11/6eda2dc2-9138-4dbc-8490-09c5c2f648f4n%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Sunnyboy010101]

On 3/12/2022 1:02 PM, Rob Pike wrote:
> Indeed. Be careful about the code Numerical Recipes if you are doing
> important numerical work. My numerical colleagues warned me off it
> years ago, as the algorithms there do not worry enough about accuracy
> and stability. It would not be described by them as a "great book".
>
> For simple stuff it's fine, of course, but floating point is a
> minefield that requires expert navigation.

Indeed. I was writing reservoir simulators in the 80s in FORTRAN.  When
we moved from a Cyber CDC  68bit word machine to a Honeywell MULTICS
smaller word machine, we had to spend months retooling the numerical
sparse matrix solvers because they would not converge on the smaller
double-word MULTICS FORTRAN.

>
> -rob
>
> On Sun, Mar 13, 2022 at 4:20 AM Neal G. <cven...@earthlink.net> wrote:
>> Mark,
>>
>> Yes, [6,30]HPSORT.FTN is the algorithm I had used, revised to sort an integer array, and to replace the IF blocks with GOTO so I can compile with older versions of FORTRAN.
>>
>> Similar to your example, my test program measures the sort time as a crude performance gauge between different systems, and a comparison with an alternate sort algorithm I had used many years ago.
>>
>> I swapped the HPSORT.FTN code into my test program and confirmed the problem is present in their code and unrelated to my revisions.
>>
>> The problem is an unstated limitation, a deficiency, in that Numerical Recipes algorithm. The maximum size of the array the algorithm can sort is 1/2 the maximum value of INTEGER. If this limit is exceeded, one of the index variables will overflow and appear to be negative on the next pass through the algorithm. The algorithm does not detect this out of bounds value and attempts to use the negative index to access the array.
>>
>> In my test program the array size was 20000. Larger than 1/2 the max value of INTEGER*2. The default DEC F77 INTEGER. So when run, this problem occurs and the memory mapping mechanism which implements the DEC F77 Virtual Array caught this invalid memory access and terminated the program.
>>
>> When compiled with /CK (DCL /CHECK) the out of bounds situation is detected and reported clearly.
>>
>> Oddly, when I had used the test program with normal arrays, no error was detected; even which compiled with /CK.
>>
>> Curious, if the sort was actually working I put the sorted values into a spreadsheet and plotted the result. Yes, even though the sort appeared to complete without error, about one quarter of the values were not sorted correctly.
>>
>> I'm not entirely certain of the moral of the story; be wary of example code, maybe. But it is certainly enlightening to dive into some of the basic issues developers faced on these classic systems.
>>
>> - Neal G.
>>
>> --
>> You received this message because you are subscribed to the Google Groups "[PiDP-11]" group.
>> To unsubscribe from this group and stop receiving emails from it, send an email to pidp-11+u...@googlegroups.com.
>> To view this discussion on the web visit https://groups.google.com/d/msgid/pidp-11/6eda2dc2-9138-4dbc-8490-09c5c2f648f4n%40googlegroups.com.

--
This email has been checked for viruses by Avast antivirus software.
https://www.avast.com/antivirus

[Neal G.]

Johnny,

Here's the source, using the HPSORT routine; and what I see when run using a normal array and a virtual array.

$ TYPE UHEAP4X.FTN
      PROGRAM UHEAP4

      REAL Q,T1,T2
      INTEGER II,ZDATA,SZ
      INTEGER*4 JS

      DIMENSION ZDATA(20000)
      SZ=20000
      JS=719

      WRITE(5,105) SZ
  105 FORMAT(1X,11HItem count: ,I6)

      DO 125 II=1,SZ
      ZDATA(II) = INT( RAN(JS) * 10000.0 )
  125 CONTINUE

      T1 = SECNDS(0)
      CALL HPSORT(SZ,ZDATA)
      T2 = SECNDS(0)

      T2=T2-T1

      WRITE(5,140) T2
  140 FORMAT(1X,'Sort time: ',F12.6)

      END

C ---------------------------------------------------------
C -- Numerical Recipes heap sort example revised to
C -- sort an integer array
C --
      SUBROUTINE hpsort(n,ra)

      INTEGER n,ra,i,ir,j,rra
      DIMENSION ra(n)

      if (n.lt.2) return
      l=n/2+1
      ir=n

410   continue

        if(l.gt.1)then
          l=l-1
          rra=ra(l)
        else
          rra=ra(ir)
          ra(ir)=ra(1)
          ir=ir-1
          if(ir.eq.1)then
            ra(1)=rra
            return
          endif
        endif
        i=l
        j=l+l

420     if(j.le.ir)then

          if(j.lt.ir)then
            if(ra(j).lt.ra(j+1))j=j+1
          endif
          if(rra.lt.ra(j))then
            ra(i)=ra(j)
            i=j
            j=j+j
          else
            j=ir+1
          endif

        goto 420
        endif
        ra(i)=rra
      goto 410
      END
$
$ FORTRAN/Check UHEAP4X
$ LINK UHEAP4X,LB:[1,1]F77FCS.OLB/LIB
$ RUN UHEAP4X
Item count: 20000
Sort time:    40.937500
$

... edit uheap4x.ftn to replace DIMENSION with VIRTUAL ...

$ FORTRAN/Check UHEAP4X
$ LINK UHEAP4X,LB:[1,1]F77FCS.OLB/LIB
$ RUN UHEAP4X
Item count: 20000

TT3  --  ERROR 94
Array reference outside array

  in   "HPSORT" at or after 22

  from "UHEAP4" at or after 13

TT3  --  ERROR 112
Virtual array mapping error

  in   "HPSORT" at or after 22

  from "UHEAP4" at or after 13

TT3  --  Exiting due to ERROR 3
Odd address trap (SST0)

at PC = 003076
  in   "HPSORT" at or after 22

  from "UHEAP4" at or after 13

$

- Neal G.

[Johnny Billquist]

Thanks. Yes, it definitely seems /CK isn't doing its job.

The second time the code hits the line:
if(ra(j).lt.ra(j+1))j=j+1
in hpsort, j is -25540, which is clearly out of bounds. Caught if
virtual, but if dimension...


Hmm. I think I see where/why this fails to be detected. The problem is
that F77 isn't checking an out of bound access by checking indexes as
such. Instead it computes the offset based on the indexes and data size,
and then checks if this address computation leads to a reference outside
of the memory space of the array referred to.

So, basically, in this case:

RA is an array with 20000 elements (of integer*2), which means it has an
address allocation of 40000 bytes.

The index -25540 is multiplied by the size of elements (integer*2),
leading to 14656 as the offset (it's all calculated as words, and
unsigned at that). -25540 is 40096 unsigned. Doubled, and truncated to
16 bits is 14656.
Now, 14656 is within the 40000. /CK thinks all is good.

With the virtual array, however, addresses are in fact calculated with
more than 16 bits, since now the limit is just that there can be at max
65536 elements in any array, independent of data type. And it's mapped
into a memory region, which can be pretty big. So, again, -25540 is
(unsigned) 40096. Multiplied by 2 is 80192. As we are not truncating to
16 bits here, this is clearly outside the 40000, and you get the error.

This was interesting. Thanks for bringing it up. I learned something I
hadn't thought about.

Johnny

> --
> You received this message because you are subscribed to the Google
> Groups "[PiDP-11]" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to pidp-11+u...@googlegroups.com
> <mailto:pidp-11+u...@googlegroups.com>.
> To view this discussion on the web visit

> https://groups.google.com/d/msgid/pidp-11/68142db0-d732-4a63-83e1-a31f7b281463n%40googlegroups.com
> <https://groups.google.com/d/msgid/pidp-11/68142db0-d732-4a63-83e1-a31f7b281463n%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Mark Matlock]

Johnny,
Thanks for the explanation! So maybe the moral of the story is that if you are working on arrays that have indexes that you think might be causing a problem you should make the arrays virtual to take advantage of more careful index checking by F77’s /CK switch?

Also, in my experience the PDP-11 symbolic debugger has helped me find obscure errors when I couldn’t figure it out otherwise.

Best,
Mark

> To unsubscribe from this group and stop receiving emails from it, send an email to pidp-11+u...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/pidp-11/61bce367-0347-5976-57eb-857c108be920%40softjar.se.

[Johnny Billquist]

Using virtual do allow you to catch some out of bound errors that F77
won't catch with usual arrays, yes. However, it mainly becomes a thing
if you just have random crazy index values, and are a bit unfortunate.
But the larger the array is, the higher is the risk that you are unlucky
and /CK isn't catching the problem, in which case virtual can help.

If you are just running through indexes, then you'll be caught as soon
as you step out.

By the way, because of the way the checking works, be even more aware of
the problem of multidimensional arrays. And virtual won't help there.
The non-dominant index in an array can be seriously out of range without
being caught.

The symbolic debugger is a different thing/tool, which is a very good
help as well. I just wish I had some sources to it, so it could be
expanded to more languages. It only supports Fortran and Cobol, and of
course just doing things in pure Macro-11, which works for any language.
But it would be cool if it could be taught other high level languages as
well... My only other complaint about the symbolic debugger is about
performance, and the slight problem of debugging programs that do
advanced screen handling.

Johnny