Challenging Function

[JohnNichols]

FUNCTION   TABA (NPC,SF,T,PA)                                     TABS 915
      DIMENSION PA(2,NPC)                                               TABS 916
C                                                                       TABS 917
C     SPECTRUM INTERPOLATION                                            TABS 918
C                                                                       TABS 919
      DO 100 I=2,NPC                                                    TABS 920
      T1=PA(1,I-1)                                                      TABS 921
      T2=PA(1,I  )                                                      TABS 922
      IF(T.LE.T2) GO TO 200                                             TABS 923
  100 CONTINUE                                                          TABS 924
  200 R1=(T2-T)/(T2-T1)                                                 TABS 925
      R2=(T-T1)/(T2-T1)                                                 TABS 926
      TABA=SF*(PA(2,I-1)*R1+PA(2,I)*R2)                                 TABS 927
      RETURN                                                            TABS 928
      END                             

In modern Fortran I is equal to five as NPC is 4, any one got any ideas what value I would have had in Fortran 66 at the end of a do loop?

I will need to look at the original formula as well.

[DavidWhite]

There are some comments in here:

http://www.ibiblio.org/pub/languages/fortran/ch2-18.html

[LRaim]

Very probably the value of I was I when the DO was terminated by the GOTO 200.

When exiting through 100 CONTINUE the value could have been either I=NPC or I=NPC+1 depending on the compiler implementation of the DO loop.

The code can be considered  CORRECT if T is always LE PA(1,NPC) or WRONG if this constraint is not verified. But the programmer could have verified the value of T before calling the TABA subroutine.

[TimP]

This fragment would have been non-portable in Fortran 66, leaving I undefined.  The result might even change with compiler flags in the same compiler.  The last f66 compiler I used had 3 different possibilities with different compile flags: 0, NPC, NPC+1.  Simply adding I=NPC+1 after 100 CONTINUE would have made it work the same as it would after standardization.

The bigger problem is of course that T1 and T2 are undefined for the case NPC==1, although a majority of f66 compilers would have set them to PA(1,1) and PA(1,2).

Even in Fortran 77 with the KAP pre-processor the result would differ according to compile options.

It was interesting to see David's quotation of my post from the time before I worked in the computing industry.

[mecej4]

My suggestion is that, rather than trying to reproduce the quirks of, say, the Fortran compiler on an IBM 704, you should replace the function by one in modern Fortran that performs the intended function: linear interpolation in a table.

The first row of array PA should be sorted in ascending order and, unless extrapolation is acceptable, the interpolation point T should not lie outside the interval [PA(1,1), PA(1,NPC)].

[TimP]

We had a 704 still in production when I was first employed after completing my physics degree, and you could have bought replacement vacuum tubes with my employer's logo on them.

[LRaim]

But the function TABA does perform linear interpolation.
​To be safe, one should only add "I1 = I" in the do loop and use I1 outside the do.

[Les_Neilson]

Luigi R. wrote:

But the function TABA does perform linear interpolation.
​To be safe, one should only add "I1 = I" in the do loop and use I1 outside the do.

Actually if you are using the code as shown and expect I to have the value NPC outside of the do loop, you could just use NPC and NPC - 1 instead.

<edit> Oh forget that. I missed the goto 200 which means you would not necessarily want the values at NPC and NPC - 1 Sorry. Yes use an intermediate variable I1 instead. </edit>

Les

[JohnNichols]

mecej4 wrote:

My suggestion is that, rather than trying to reproduce the quirks of, say, the Fortran compiler on an IBM 704, you should replace the function by one in modern Fortran that performs the intended function: linear interpolation in a table.

The first row of array PA should be sorted in ascending order and, unless extrapolation is acceptable, the interpolation point T should not lie outside the interval [PA(1,1), PA(1,NPC)].

I agree it must be recoded -- I found this yesterday and was intrigued by the result and it was late and I was tired of hunting these errors. Clearly it needs to be rewritten and that is not difficult - although in reality  I am not a big fan of the spectrum method - it was essentially invented in 1916 and is still used a lot. We have much better methods.

I was perhaps using this as example of some of the interesting coding issues in this program compared to some of the others from Powell's group.  It is the topological ideas that are advanced here, the code is not sophisticated at all.  (https://en.wikipedia.org/wiki/Edward_L._Wilson​) is the author I believe, but he was not a grad student at the time but on faculty most likely. It might explain his stumbling Fortran compared to the grad students who appeared to share code and ideas.

He is a famous engineer.

[JohnNichols]

Actually there are bigger issues coming into this function, T has a strange set of values coming into the function which I have to solve first, including some negative sqrt's that need to be trapped.

[JohnNichols]

The IBM 704 had a 38-bit accumulator, a 36-bit multiplier quotient register, and three 15-bit index registers. The contents of the index registers were subtracted from the base address, so the index registers were also called "decrement registers". All three index registers could participate in an instruction: the three-bit tag field in the instruction was a bit map specifying which of the registers would participate in the operation. However, when more than one index register was selected, then their contents were or'ed – not added – together before the decrement took place. This behavior persisted in later Scientific Architecture machines (such as the IBM 709 and IBM 7090) until the IBM 7094. The IBM 7094, introduced in 1962, increased the number of index registers to seven and only selected one at a time; the "or" behavior remained available in a compatibility mode of the IBM 7094.^[11]

[JohnNichols]

Interesting the Wikipedia site says FORTRAN and LISP were invented for the 704.

I have been using the MKL Random Number Gaussian Routines also as part of this larger research project on bridge vibration.  The gold standard is the Box Muller methods from Princeton circa 1958 - it is in MKL.  Our librarian found one of the minor but not complete references from the INTEL routine documents, but we could not find Technical Report 13 anywhere which is the detailed report for the Army that forms the basis for the MKL routine.  We spent a few days talking to the Princeton Library and somewhat the math department. Finally we were sent No 9, which has been scanned and is attached.  Not the real one - but close.

We were asked by the Princeton Librarian in a round about sort of way - should we scan Technical Report 13.  TR13 is not available anywhere at the moment hopefully we can get it stored safely on the web, which is a lot safer than a book in a library (weird )

I replied as loudly as possible in an email YES.

The reports discuss quite clearly the impact of the 704 on this mathematics and the uses to which it was put.

Interestingly the report author thanks a Tukey for help and Mrs A Schay for the math.

There are more computers in my building than the total number of 704's that were ever sold.

My father in law just told me Purdue did not have a 704

--------------------------------------------------------------------------------------------------

No, in that era Purdue had a Univac SS80, then an IBM 7094.

I heard a lot about the IBM 704, but never used one.

Where did you come across that memory?

Vic

[John_Campbell]

BAD DATA !

This is an interpolation routine from TABS or ETABS
If T (time) > PA(1,NPC), in the original pre F77 version "I" would have the value NPC rather than NPC+1
What this situation means is that the value T is out of range for the data values supplied in PA(1,1:NPC) so you should check the data you have provided. Isn't this data for a time varying force or acceleration and your time range exceeds the defined data.

The supplied data should typically be:
PA(1,1) = 0
PA(1,NPC) = large time value
PA(2,1) = 0
PA(2,NPC) = 0

Some versions have error tests for the data to be in range.

[JohnNichols]

Thanks

Yes - I tracked down the error and put in checks.  It now returns the correct value.