I'm writing a toy Fortran 66 interpreter for fun -- want accuracy

Help
41

The following are a Fortran IV Gauss-Jordan elimination routine that computes the inverse of a matrix along with the magnitude of its determinate and the det mag. scaled by the Euclidean norm of the matrix which I guess is a measure of how well the matrix is conditioned and my modern Fortran refactoring. The source is a verbatim copy (in all its statement label and GO TO glory) of a routine listed in Hornbeck’s “Numerical Methods” book from 1975.

I tested with ifx 2025.3 and gfortran-13 and both give the same results for the following matrix

  1.0  3.0  3.0  2.0  1.0
  1.0  4.0  3.0  3.0 -1.0
  1.0  3.0  4.0  1.0  1.0
  1.0  1.0  1.0  1.0 -1.0
  1.0 -2.0 -1.0  2.0  2.0

which has as its inverse (to 1 decimal place output)

  2.0 -1.3 -1.0  1.7 -0.3
  2.0 -0.7 -1.0  0.5 -0.3
 -2.0  0.8  1.4 -0.6  0.4
 -1.0  0.8  0.4 -0.6  0.4
  1.0 -0.5 -0.4 -0.1 -0.1

   det mag = 15.0

This is the Fortran IV code

      SUBROUTINE INVDET(C,N,DTNRM,DETM)
      DIMENSION C(70,70),J(120)
      PD=1.
      DO 124 L=1,N
      DD=0.
      DO 123 K=1,N
  123 DD=DD+C(L,K)*C(L,K)
      DD=SQRT(DD)
  124 PD=PD*DD
      DETM=1.
      DO 125 L=1,N
  125 J(L+20)=L
      DO 144 L=1,N
      CC=0.
      M=L
      DO 135 K=L,N
      IF ((ABS(CC)-ABS(C(L,K))).GE.0.) GO TO 135
  126 M=K
      CC=C(L,K)
  135 CONTINUE
  127 IF(L.EQ.M) GO TO 138
  128 K=J(M+20)
     J(M+20)=J(L+20)
      J(L+20)=K
      DO 137 K=1,N
      S=C(K,L)
      C(K,L)=C(K,M)
  137 C(K,M)=S
  138 C(L,L)=1.
      DETM=DETM*CC
      DO 139 M=1,N
  139 C(L,M)=C(L,M)/CC
      DO 142 M=1,N
      IF(L.EQ.M) GO TO 142
  129 CC=C(M,L)
      IF (CC.EQ.0.) GO TO 142
  130 C(M,L)=0.
      DO 141 K=1,N
  141 C(M,K)=C(M,K)-CC*C(L,K)
  142 CONTINUE
  144 CONTINUE
      DO 143 L=1,N
      IF (J(L+20).EQ.L) GO TO 143
  131 M=L
  132 M=M+1
      IF (J(M+20).EQ.L) GO TO 133
  136 IF (N.GT.M) GO TO 132
  133 J(M+20)=J(L+20)
      DO 163 K=1,N
      CC=C(L,K)
      C(L,K)=C(M,K)
  163 C(M,K)=CC
      J(L+20)=L
  143 CONTINUE
      DETM=ABS(DETM)
      DTNRM=DETM/PD
      RETURN
      END

Here is my refactored code (much longer because I added whitespace to improve readability)

invdetmf.f90 (1.9 KB)

and a test program. Fortran IV must be compiled separately due to fixed field format.

testinvdet.f90 (3.6 KB)

Also a question for the mods. Why can’t I upload a file with a .f file extension. With my browser (firefox) it only appears to support .f90

Edit. For some reason the original code required the J array dimension which is used for column pivoting to be at least N+21. I have no idea where this comes from.

42

Thanks @m_b_metcalf and @rwmsu!

@Beliavsky yes, it can be revived using LLMs, but I think more robust is to actually compile it, and then let the compiler produce a modern version, which I want to eventually do also with LFortran.

@themos great find! It’s interesting that they list two versions of the program, one in “English units”, one in “metric units”. Also, in their printout, it seems I need to use the following substitutions: # -> =, and % -> ( and ◻ -> ). I still need to figure out what Ɛ means. Does anyone know why they would use such a character set?

43

@m_b_metcalf I got your program to work using GFortran by making two sets of changes:

  • READ INPUT TAPE 5, 51, ANAME, N → READ(5,51) ANAME, N (it looks like read input tape was introduced in FORTRAN II, but changed in FORTRAN IV to the modern read() syntax)
  • SQRTF → SQRT

Then it compiles and runs:

$ gfortran -std=legacy a.f
$ cat input
TEST  1
      1.00     -5.00      6.00
$ ./a.out < input
1ROOTS OF QUADRATIC EQUATIONS FROM  ST
0SET NO.  1
 A =     1.00            B =    -5.00            C =     6.00
 REAL                          X1 =      3.000             X2 =      2.000
0END OF   ST

The above solves the equation x^2 - 5x + 6 = 0 with roots x=2 and x=3.

It looks like the A Fortran Primer book is from 1963.

Does anyone know if any compiler supports the read input tape and write output tape syntax? I don’t know if it is worth it, but if it is just a few legacy extensions like this, then it would be very nice to be able to compile these old codes unmodified.

45

Please don’t. :grinning_face_with_smiling_eyes: Then you’d have to implement the FORTRAN II version of COMMON blocks and EQUIVALENCE - which worked somewhat differently than FORTRAN IV and Fortran 66.

If you look at the IBM Fortran IV manuals over at bitsavers, there were a number of things that did not go into the Fortran 66 standard. Most did go into F77 - in one form or another. So even though folks may consider FORTRAN IV and Fortran 66 the same, they aren’t.

The Bell Labs folks claim they had the first full f77 compiler.

The CDC FTN5 compiler was not far behind them. The FTN5 Reference Manual revision history says the original release was 7/20/1979. The first F77 compiler I used was the U of Minnesota M77 compiler for CDC machines in probably 1980 or 1981. M77 was an extension of their MNF compiler.

Across the river from CDC at Cray, the CFT compiler for the CRAY-1 and X-MP gradually added F77 features and was feature complete by 1983 or so. On the CRAY-2 side, there was a fork prior to this. So the initial version of CFT2 didn’t have CHARACTER data type. The Cray CFT77 compiler was, of course, feature complete on all systems. CFT77 was initially written in PASCAL - as there wasn’t a C compiler for Crays when the project started.

IIRC, IBM took a long time to release a Fortran 77 compiler.

Sperry Univac introduced their “ASCII” Fortran 77 compiler for their 36-bit computers in the early '80s. As part of it they moved from their 6-bit FIELDATA character set to 9-bit ASCII. This caused a fairly major amount of work by their customers. Interestingly, a bunch of the Cray CFT77 compiler group were refugees from Univac.

46

The ORGANICK book was published in 1963-1965. I would say that a decade later, by 1975 or so, most textbooks would instead present the stable way to compute the two real roots. Using modern fortran syntax, this is

REAL1 = -(B + SIGN(SQRT(D),B)) / (2.0*A)
REAL2 = C / (A * REAL1)

I do not know how this would have been written in FORTRAN II. The key point in the first statement is to always add terms that have the same sign; that eliminates unnecessary roundoff that would occur when adding two terms with different signs. Then the other root can be computed with the division. A second SQRT() is not required. In a more advanced discussion, special cases would be considered for small values of A or large values of B, which lead to the use of additional formulas to avoid unnecessary overflow.

47

It must be a ‘+’.

48
Substitution Original Confused with
& + *
# = -
% ( )
□ ) (

If it’s deliberate, then one guess is that = and - can be mistaken for each other, and the same goes for + and *.

49

Character sets were a disaster back then. They were a true example of “the nice thing about standards is there are so many to choose from”. So much legacy and re-invention of the wheel, from Herman Hollerith, to 5-level Teletypes, to a plethora of 6-bit codes, to EBCDIC and ASCII. ASCII finally rationalized things. (EBCDIC is weird because IBM just added a couple leading bits to one of their BCD character sets. This led to discontinuities in the A-Z and a-z encodings.)

It is why old reference manuals almost always have an “Appendix A” which shows the character encodings the compiler or system supports.

This page shows some of the sad state of affairs:

50

The continued responses here are great to see, and have me thinking about extending my interpreter with more functionality and for more use-cases once my current project is done.

One thing I notice about the Fortran 66 spec is a 6-character identifier (and 5 for Basic Standard Fortran). From what I’ve seen, implementations would typically allow for more characters to be present in ids, they just wouldn’t be used to make the id unique. Is that true? The compiler would consume any extra letters/digits, but then toss them away and just use the first 6? Or were they more strict and give errors for ids longer than 6? My practical reference here is that gfortran using –std=legacy accepts ids longer than 6 characters long (and uses all characters for uniqueness).

51

From all the responses here, it seems it makes sense to get a good compiler (and interpreter, thanks @csmoak!) going for F66, but going earlier is problematic because codes used various extensions that were different for each compiler. That being said, I do want to run some of these old codes, and not via AI “translation”, but actually using a compiler. My motivation is to see if it is possible and how backwards compatible Fortran really is, and also for marketing of Fortran I want to know if with minimal modifications we can run these old codes today. There is academic value in this “archeology” too. Some of these reports are examples of first “reproducible papers”, they give you self-contained code that you can run (well, we’ll see!).

52

Probably back then you would manually convert to punch cards anyway, so it would be easy to do this translation by hand. You can see that the quality of the pdf document is not super great, so it’s actually helpful that they take care of some of the printing quality issues this way.

53

As already noted by others, this was due to multiple incompatible character set mappings. In my first programming course and first numerical methods course that I took in school (1973 and 1974), the students had access to an old card punch. When we wrote fortran code on that card punch we had to do those types of substitutions. I remember in particular using % for a parentheses (I don’t remember now if it was left or right). When those cards were read by the compiler and printed on the line printer they came out right. There were newer card punches where this was not required, but they were reserved for the college staff and faculty. In the next semester, that old card punch was replaced and the newer ones after that did not require us to make those character substitutions.

Years later when I read the fortran standard, I noted that it specifically excluded how the native character set on a processor was translated to the fortran character set. The f77 text for the list of exclusions is in section 1.3.2. This includes, “The method of transcription of programs or their input or output data to or from a data processing medium.” Because of my experience with that card punch, I knew what that meant and why it was in that list.

54

@RonShepard I see, that makes perfect sense. Also probably not a coincidence that you also had the same substitution % → (.

55

Limited Character Sets (CDC/Mainframe Era): Early Fortran (particularly on Control Data Corporation - CDC - machines) used a restricted 64-character set to ensure compatibility across printers and punch card systems. In this set, symbols often had to serve multiple purposes, and the equals sign was sometimes replaced by the # symbol for assignment statements. ASCIIi 6/12 and Bang!Bang! were particularly interesting looking back. I wrote something called “EXODUS” that migrated hundreds of users from older machines to our first ASCII machines and did a LOT of file conversions automatically; expanding UPDATE libraries back to source, converting ASCII 6/12 to ASCII, routines to read COS and SCOPE binary files, dealing with fixed width files versus variable length files, line terminators, multi-file files, … Not sure about your entire list but at least two of those ring a bell and make me think an older CDC machine was involved.

Punch Card Limitations: Before modern keyboard standards, characters on a card were punched as a series of holes. If too many holes were punched close together, the card might tear. The character set was limited to prevent this, which caused, for example, certain special characters to be substituted.

56

ALGOL faced similar design issues related to translation between mathematical language and the hardware representation. As a result they defined three languages:

  • Reference Language
    1. It is the working language of the committee.
    2. It is the defining language.
    3. The characters are determined by ease of mutual understanding and not by any computer limitations, coders notation, or pure mathematical notation.
    4. It is the basic reference and guide for compiler builders.
    5. It is the guide for all hardware representations.
    6. It is the guide for transliterating from publication language to any locally appropriate hardware representations.
    7. The main publications of the Algol language itself will use the reference representation.
  • Publication Language
    1. The publication language admits variations of the reference language according to usage of printing and handwriting (e.g. subscripts, spaces, exponents, Greek letters).
    2. It is used for stating and communicating process.
    3. The characters used may be different in different countries, but univocal correspondence with reference representation must be secured.
  • Hardware Representations
    1. Each of these is a condensation of the reference language enforced by the limited number of characters on the standard input equipment.
    2. Each one of these uses the character set of a particular computer and is the language accepted by a translator for that computer.
    3. Each of these must by accompanied by a special set of rules for transliterating from publication or reference language.

More about the hardware representations can be read in,

GCC16 will have a new ALGOL 68 frontend: Top (The GNU Algol 68 Compiler)

57

Do they really? News to me. I found a web page claiming that Bell Labs developed Fortran, PL/I and System\360 Fortran - clearly nonsense, and from the initial part of the page it seems it is AI-generated slop. To the best of my knowledge, Bell Labs didn’t have a Fortran compiler at all..They did create f2c, which is not a Fortran compiler, and was developed many years after F77 compilers went mainstream.

If you have evidence otherwise, I’d love to see it. I’m not sure if DEC was first with a full FORTRAN 77, but if not, it was second.

58

August 1, 1978:

59

Thanks - I had never heard of that and didn’t run across anyone using it.

60

CDC released theirs in 1979. So if VMS 2.0 was fully F77 compliant in 1980, DEC was at best third. IBM VS Fortran came out in Feb 1981. Not sure when the Univac ASCII compiler came out. The manual at bitsavers is from 1982, but there is no revision history - so might not be the first release.

61

I knew about this compiler, but only later, maybe 1982 or so. My division had a VAX 11/780, and we considered whether we should run the VMS or the unix operating system on it. One piece of evidence used in that decision was to compare runtimes for some benchmark programs with the DEC f77 compiler on VMS or the Bell Labs f77 compiler on unix. The DEC compiler was much faster, maybe a factor of 3x or 4x, I don’t remember exactly, so we decided to use the DEC software on that machine. I think the intention at that time was to revisit that decision periodically, but we ended up using DEC software on that machine until it was replaced (which I think was maybe 1994 or so). Another division did run unix on their VAX, so we could experiment with codes and compilers on that other machine whenever we wanted.

That Bell Labs f77 document is interesting. It describes the low level implementation of several f77 features, such as the hidden dummy arguments for character variables and the use of _ for fortran external symbols. It uses implicit undefined, rather than what I think was the more common extension implicit none. I knew a few other compilers that had that extension, and I suspect that they copied the idea from the Ball Labs compiler. It also allowed both " and ' to delimit character strings, with what looks like the later f90 conventions. I don’t remember any other f77 compilers that did that until much later, after the fortran 8x drafts had been circulating. I have also read elsewhere that the f2c cross compiler was based on this Bell Labs f77 compiler conventions.