Community
cancel
Showing results for
Did you mean:
Highlighted
New Contributor II
79 Views

## Challenging Function

```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.

Tags (1)
13 Replies
Highlighted
Black Belt
79 Views

There are some comments in here:

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

Highlighted
New Contributor I
79 Views

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.

Highlighted
Black Belt
79 Views

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.

Highlighted
Black Belt
79 Views

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)].

Highlighted
Black Belt
79 Views

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.

Highlighted
New Contributor I
79 Views

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.

Highlighted
Valued Contributor II
79 Views

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

Highlighted
New Contributor II
79 Views

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.

Highlighted
New Contributor II
79 Views

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.

Highlighted
New Contributor II
79 Views

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]

Highlighted
New Contributor II
79 Views

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

Highlighted
New Contributor II
79 Views

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.

Highlighted
New Contributor II
79 Views

Thanks

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