topic Re: Harrison 1967 Fortran in Intel® Fortran Compiler

Harrison 1967 Fortran

JohnNichols — Fri, 16 Jan 2026 15:02:46 GMT

Do 610 I=1,L IP1 = I + 1 TEMP = ABS(ASAT(I,I)) K = I Do 520 J = I,L DUM = ABS(ASAT(J,I)) - TEMP IF(DUM)520,520,510 510 K = J TEMP = ABS(ASAT(J,I)) 520 Continue IF(K - I) 530, 550,530 530 DO 540 J = I,KJ TEMP = ASAT(I,J) ASAT(I,J) = ASAT(K,J) ASAT(K,J) = TEMP 540 END Do 550 IF(ASAT(I,I)) 570, 1000, 570 570 TEMP = 1.0/ASAT(I,I) do 580 J = I,KJ ASAT(I,J) = ASAT(I,J)*TEMP 580 END DO do 610 J = 1,L IF(I-J)590,610,590 590 TEMP = ASAT(J,I) DO 600 K = IP1,KJ ASAT(J,K) = ASAT(J,K) -TEMP*ASAT(I,K) 600 END DO 610 END DO DO 650 I = 1,M3 CSAT(I) = 0.0 DO 651 J = 1,L CSAT(I) = CSAT(I) + A(J,I)*ASAT(J,L+1) 651 END DO 650 END DO! DO 660 I = 1,M2 K = ((I+1)/2)*2-1 SATX(I) = SF(I,1)*CSAT(K) + SF(I,2)*CSAT(K + 1) 660 END DO

This comes from some Harrison Fortran Code in 1967 for second order elastic analysis, he says that it ran on a CD 3600 computer in 17 minutes for 17 cycles, the code write-up does not really discuss this section, after playing with it for a while, I think this is the earliest inverter I have ever found.

What was a CD 3600 computer? It was at Le High University and would 17 minutes be the limit they were allowed for time. I put in a separate inverter just to make sure from Conte and De Boore, as I tried to work out the code. They really worked to use limited memory.

Why is this important, running on modern computers it allows us to do major structural problems not constrained by the annoyances in the major programs. I was reading Strand7 stuff on this type of program and it was so close to some of the stuff Harrison wrote it was eerie.

It runs in a few seconds.

Editing out the arithmetic ifs is fun.

Re: Harrison 1967 Fortran

mecej4O — Fri, 16 Jan 2026 15:36:11 GMT

Control Data Corporation sold 36xx model computers prior to creating its supercomputer line, the 6600/6400 models, in the 1960s. The CDC 6000 was in many senses the world's first supercomputer.

Re: Harrison 1967 Fortran

JohnNichols — Sat, 17 Jan 2026 16:01:53 GMT

Thank you. Now I remember using a CDC machine at a computer place in the early 1980's before we got PC's. I used to run ULARC and AXISHELL to design steel structures for coal mines. These programs came from Powell's group at UCB. Brings back memories, we had the programs on a tape that my boss bought at UCB when he was a grad student. It included TABS, but I had a a lot of trouble getting that to run and it was not really useful.

Back to Fortran, Harrison and a few others developed elastic and plastic structures programs in the late 1960's and early 1970s. His book and papers are really helpful to understand why they did what they did. Essentially without great data they designed elastically and the ventured into plastic using Fortran. It appears that the first inverter was written by Wang in the mid 1960's I am trying to get a copy of the paper.

But the current problem is we now have great longitudinal frequency data from steel and concrete bridges, and the results mirror the interesting points that Harrison and the others raised about the problems they saw with the analysis, but lacked the data to analyze. Plus they were just getting the programs accepted.

As an aside it is sad to see civil engineering depts teaching this stuff with MatLab, you cannot design a real structure with MatLab.

The question is can we amend a Fortran program based on plastic to include strain hardening, the need is the elastic/plastic analysis on old structures has a probability of failure approaching the lower limit, the frequency data leads to this refined result, but the structures are still "ok". It has to be the strain hardening, I think.

Can I use Newton Ralphson to amend the plastic program to follow the strain hardening increase? If I use Harrison's second order analysis I can get the worst case strain e and then look for the capacity, which is greater than the plastic capacity. Or is this a waste of time?

Re: Harrison 1967 Fortran

mecej4O — Sat, 17 Jan 2026 16:58:23 GMT

From what I read from the page extracts of the Harrison book that you posted, what they did in the book is to account for geometric nonlinearities, not plastic stress-strain relationships. I have never used commercial FEM software, but I recall that my colleagues who worked in solid mechanics considered ABAQUS to be the best software for modelling nonlinear elasto-plastic problems. There are a number of people who post in this forum who are users of ABAQUS, let's hope that they add to this thread.

Here is a link to photograph of a CDC computer being unloaded from an airplane:

The CDC Computer System Reaching Bombay — Google Arts & Culture

Re: Harrison 1967 Fortran

JohnNichols — Sat, 17 Jan 2026 21:52:48 GMT

https://caeassistant.com/product/abaqus-kinematic-hardening-1/#1674478660080-7f8862a6-aa43

Yes, it is good software, as is the European Dlubal software and Strand 7, and I have used all of them, but they do not fit inside a Monte Carlo simulation to do 2,000,000 iterations in a night for a Structural Reliability Problem. And getting the answers out can be a pain, even with Fortran. At a constant thermal strain, EI and EA are not constant with time, but linearly decreasing. Consider the structural eigenvalue equation with EI dropping with time, then the frequency drops and we are now measuring that with t-Stats in the hundreds.

But, their outline provides a path forward for this type of modelling, thanks.

From Abaqus - this is one model, which is the alternative to a simple two element linear model.

Re: Harrison 1967 Fortran

JohnNichols — Sat, 17 Jan 2026 22:19:37 GMT

Here is a paper I just got, the opening paragraphs show the problem.

You have one day to determine the structural reliability beta for a bridge, you are not going to get there with Abaqus, if you could there would be papers on this matter and the people I work with would be discussing it.

The main structural frequency is dropping at 1500 microHz per day with a t-Stat of say 50 on the linear regression elements, we are trying to pattern match these curves. The best solution is Fortran with a special purpose program.

Re: Harrison 1967 Fortran

witwald — Sun, 18 Jan 2026 03:21:16 GMT

@JohnNichols wrote:
This comes from some Harrison Fortran Code in 1967 for second order elastic analysis, he says that it ran on a CD 3600 computer in 17 minutes for 17 cycles, the code write-up does not really discuss this section, after playing with it for a while, I think this is the earliest inverter I have ever found.

@JohnNichols Would you be able to share the code for the PROGRAM ELSO that was located in Chapter 11 in Harrison's textbook from 1973, "Computer methods in structural analysis"? I'd like to try and run it. The original code would be ideal, if that's what you started with.

Re: Harrison 1967 Fortran

JohnNichols — Sun, 18 Jan 2026 14:20:37 GMT

Here is the Harrison code, it has been cleaned up so it runs on easily to check the internals, i.e. implicit none, etc.. It runs in VS and the latest One API.

The input is a.inp, I added one line to make the code simpler on the arithmetic ifs.

I also added a inverter as I was in a hurry and I was having trouble with Harrison's. inverter, it kept giving me zeros on the first entry, it took a while to find the mistake. I am using model 11.9 for the run, Harrison gives an answer, and mine are similar except for the starting deflection on y at node 3, which linear gives -5 and in the end converges to -3, which is about right. Harrison says it is -2.9 linear, I have not run his other models. There are some typos in his code, I think I got all of them.

HIs original code will not run now as the input statements are for the CDC 3600.

Re the other code, it will take me a day or so to package it up.

Re: Harrison 1967 Fortran

JohnNichols — Sun, 18 Jan 2026 14:36:30 GMT

RE Abaqus, I would love to use ABAQUS for modelling, I used it for a long time on a supercomputer. Model size is limited by the CAE, time is limited by the supercomputer people, and the runs have to be batched.

The problem is you need to write a decent Fortran routine to provide the input model in a consistent format, using a GUI on a large model is a recipe for disaster always, this is why one PhD student does one Abaqus model in three years.

Then you need to write a Fortran program to get the output into sensible numbers that you can show people.

I have this for Strand7 that I run on a PC. Strand7 on a PC now runs faster than ABAQUS in 1995 on a supercomputer, I know as I was banned from the computer at Newcastle University as STRAND7 grabbed the threads and ran overnight. The HOD called me in and said I had upset a lot of people.

So I use say Strand7 to get the natural frequencies on a bridge, it will give me 100 frequencies below 100 Hz, the mems will measure say 30 and the great Harrison code that @mecej4O helped with on 3D modelling will give a closer set to the 30. It also tells me which members generate which frequency.

Vic Abell at Purdue pulled out the last CDC 3600(??) there, they cut all the cables, and it got shipped off, it is now running somewhere and the poor people had to redo the cabling. Vic said there were a lot of cables.

He did a masters on organizing the Uni Schedule, he got it running on the CDC, but it was limited by the time he allowed. I retyped the thesis and put the code into Intel Fortran, it ran in a few seconds, about 5 years ago.

The issue we are finding is that we now have excellent acceleration data on structures, but the modelling needs to catch up and the design methods for bridges use loads from WW2.

Re: Harrison 1967 Fortran

jimdempseyatthecove — Sun, 18 Jan 2026 17:14:17 GMT

John,

I haven't looked at the code bases you mentioned, but I do have a comment.

When you describe the Ramberg-Osgood model... a, b and n are material constants. (as well as Young's modulus)...

Doesn't strain hardening affect these material properties?

IOW do these need to be variables, modified during the simulation?

Jim

Re: Harrison 1967 Fortran

JohnNichols — Sun, 18 Jan 2026 18:06:12 GMT

Harrison has developed the ELSO code so that it can do all of those changes rather easily, so adding a model like the RO should not be hard, the difficulty is the changes to the interaction model for the axial to bending load, I will need to look at that some more as it will need to fit in that section of the code. In looking at @mecej4 comments today, he as usual pointed me in the right direction.

But in terms of coding and modelling I think the better answer is a power model, rather than the Abaqus models.

strain= constant (stress/E) to power m, easy to fit and code and it is only a slight variation on the coding for the stability functions Harrison used in 2D.

The challenge is fitting all this stuff into the USA LRFD methods, I looked at one paper this morning, instead of using the standard Plastic Analysis Portal Frame, they invented a new one, leaving out the critical load and adding distributed loads. There is a small group of authors you need to follow for plastic analysis and most them are English, although Harrison in USA from Le High and Melchers is from Australia via Oxford. But there is a standard model, not using that makes your results difficult to compare.

Classic example, someone developed a bridge loading from Weigh in Motion data from a WIM in the ground, no comment on the WIM ground conditions etc., they said the large loads should be discarded. No mention of Fryba's work who had solved this forty years ago on railways.

Re: Harrison 1967 Fortran

DavidWhite — Mon, 19 Jan 2026 01:28:43 GMT

I believe a CDC 6400 was the mainframe we had at Adelaide University while I was a undergraduate in 1976. It was such an advanced machine that I believe they rented out compute time to local businesses.

As undergraduates, we had to use punch machines to punch cards and submit the decks for execution. Usually the punch machine had worn out typeewriter ribbons, so you could not read the text at the top. One got adept at reading the punch marks!

Two stories from that time - two students submitted the same deck as their homework assignment - the tutor quickly spotted the same spelling mistakes in the Fortran comments and awarded them 18/20 for the assignment, shared equally as 9/20 each!

Also some post grads using mag tape input tested code to put two copies of itself onto the input queue. Eventually the queue overloaded and the machine went down. To restart the machine, the techs used the same input tape, resulting in the same crash. Only then did they restart the machine with a new tape!

Re: Harrison 1967 Fortran

JohnNichols — Mon, 19 Jan 2026 15:53:16 GMT

We had to walk about a mile to the computer support section and drop off the punch cards and wait 24 hours. My professor did an experiment with the students he gave each one different data and had them run the analysis and then hand in the answers, it was for some research, of the class I was the only one to do the wrong data. Then I found a little box computer in the library they called daemon and it took cards and ran Fortran straight away, 'twas great.

Re: Harrison 1967 Fortran

witwald — Wed, 21 Jan 2026 11:02:57 GMT

Hello @JohnNichols,

Thanks for sharing your code. I took a look at it and couldn't see any reason why it might be giving different results from those produced by Harrison.

In any event, I've recreated Harrison's ELSO code in its entirety, and it runs and produces results that are all quite close to those that were published. That's for both the linear and the nonlinear solutions, which is comforting. Interestingly, the present code converges in 6 cycles instead of the quoted 13.

It's interesting just how close that old code was to modern free-format Fortran 90.

I also recreated Harrison's model of the 2nd frame, and the solution was also a good match against the published results.

A copy of the Fortran for the ELSO.f90 program is attached, together with the input decks that I used to analyze the two frames described in Harrison's book.

Re: Harrison 1967 Fortran

JohnNichols — Wed, 21 Jan 2026 15:48:01 GMT

Hello @witwald

Is there any chance you could put up a screen shot of the 11.9 problem with the crossed braces, complete. I am interested in the elastic and second order deflection arrays, I am sorry, I do not have time to load your program in for a week or so

Yes my convergence was about 7.

it is the -5.874 that is weirdly different to Harrison's result, but the final is much closer

If I have an error I would like to track it down. Thanks

John

Re: Harrison 1967 Fortran

JohnNichols — Wed, 21 Jan 2026 15:49:19 GMT

Jim

Why would there be such a difference between the CDC Fortran and Modern Intel 64 bit Fortran.

John

Re: Harrison 1967 Fortran

witwald — Thu, 22 Jan 2026 08:04:03 GMT

Note that a copy of the results shown below was included in the header comments for the ELSO.f90 source code that I provided.

! ELASTIC SECOND-ORDER ANALYSIS OF FRAME NO. 9 ! ! RESULTS OF FIRST-ORDER LINEAR ANALYSIS FOR LOAD SET NO. 1 ! ! JOINT X-MOVEMENT Y-MOVEMENT ROTATION ! 1 0.00000 0.00000 0.00000 ! 2 0.14407 -0.00593 0.02083 ! 3 0.13853 -2.94237 -0.00021 ! 4 0.13299 -0.01072 -0.01991 ! 5 0.00000 0.00000 0.00000 ! ! MEMBER TERMINAL APPLIED MOMENTS NEAR JOINTS AXIAL TENSION ! 1 322.7800 690.4254 1 AND 2 -10.4735 ! 2 -690.7827 -1532.4095 2 AND 3 -13.2906 ! 3 1532.4099 744.3958 3 AND 4 -13.2906 ! 4 -744.0441 -392.7302 4 AND 5 -18.9191 ! 5 0.3575 0.1763 2 AND 5 -8.8161 ! 6 -0.1789 -0.3521 1 AND 4 7.5894 ! ! RESULTS OF SECOND-ORDER ANALYSIS FOR LOAD SET NO. 1 ! ! CONVERGENCE IN 6 CYCLES ! ! JOINT X-MOVEMENT Y-MOVEMENT ROTATION ! 1 0.00000 0.00000 0.00000 ! 2 0.41568 -0.01430 0.02268 ! 3 0.33982 -3.08032 -0.00054 ! 4 0.26490 -0.01879 -0.02044 ! 5 0.00000 0.00000 0.00000 ! ! MEMBER TERMINAL APPLIED MOMENTS NEAR JOINTS AXIAL TENSION ! 1 273.2821 664.1279 1 AND 2 -14.5736 ! 2 -664.5078 -1581.6520 2 AND 3 -20.2485 ! 3 1581.6519 795.8525 3 AND 4 -20.2312 ! 4 -792.5956 -445.2905 4 AND 5 -22.9470 ! 5 0.3800 0.1827 2 AND 5 -0.0005 ! 6 -0.1627 -3.2573 1 AND 4 15.2821 ! ! ANALYSIS COMPLETED FOR FRAME NO. 9 UNDER LOAD SET NO. 1

I've lost my Fortran compiler for the time being, as my Windows 11 desktop PC crashed today and it needs a complete re-install. I'll be a little while sorting that one out!

Re: Harrison 1967 Fortran

JohnNichols — Thu, 22 Jan 2026 15:21:58 GMT

Hello @witwald

double precision, parameter :: TWO = 2.0 double precision, parameter :: FOUR = 4.0

I am glad you copied in the original code, I hope you used OCR.

Using the two programs I looked backward through the answers and found I had set the parameter FOUR to a value of 2 instead of 4.

Even the worst of us make mistakes.

Here is your file in a VS Folder.

Thanks heaps.

I have asked our library to find the last book by Harrison, could be a hard ask - but we shall see.

John

Re: Harrison 1967 Fortran

JohnNichols — Fri, 23 Jan 2026 03:04:21 GMT

Thanks.

I played with your code and my code today. The ZZ test for convergence starts as 0.005, if you take it down to 0.001 it takes 70 iterations but the answers climb back toward the elastic.

Interesting code.

Sorry about the crash, it is a b when that happens as my mother used to say.

Re: Harrison 1967 Fortran

mecej4 — Fri, 23 Jan 2026 04:01:02 GMT

Witwald, Running your program with the data file that you supplied gives the expected results (which are included as comments in the source file), but the program has a bug: It assumes that the array VL is set to zero at the beginning of the program. If this zeroing is not explicitly performed in the program, the results are undefined. Add "vl(1:18)=0.0" at the beginning of the executable part of the program.