- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

`AL(I)= sqrt(XM(I)*XM(I) + YM(I)*YM(I))`

Is there a faster way to solve this problem in fortran than the long method?

Link Copied

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

First step is to structure the code such that it vectorizes.

Second step is to see if the requirements can be configured to omit the sqrt (often this can be done).

IIIIF XM and YM are double precision but AL can afford some loss in precision, create a temporary array of single precision to receive the sum product, then in a seperate vectorize loop perform the sqrt (iow using single precision sqrt on double the number of variables per instruction).

You can (if you choose) to then improve the precision using one or two Newton-Raphson iterations.

Also, if AL is going to be used as a divisor, current CPU's can produce a fast reciprocal (at loss of precision).

Additional information: https://betterexplained.com/articles/understanding-quakes-fast-inverse-square-root/

Jim Dempsey

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

*CONCRETE CYLINDER STRENGTH AT AGE 2fl DAYS B*

In some 1972 code, which is useful, the format statements have this little gem, is this OLD FORTRAN?

* * instead of " "

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Can you paste the particular lines (some before, some after) using the ...,</>,Fortran markup.

Spaces matter.

Jim Dempsey

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

```
Do 2133 I = 1,NJOIN
FXX(I)=FXX(I) + 0.1*FXX1(I)
FYY(I)=FYY(I) + 0.1*FYY1(I)
FMM(I)=FMM(I) + 0.1*FMM1(I)
2133 CONTINUE
6013 IF(INDEXX .GT. 0) GO TO 8147
IF(NM1 .GT. 0) GO TO 4373
1975 CONTINUE
WRITE(*1976)WLIVE,WIND '
1976 FORMAT(40X,*SUPERIMPOSED LOAD=>*,F10.3,10X,*WIND LOAD=.*,F10.3)
WRITE(6,1987) \ • , \
1697 FORMAT(/// ,50X,*APPLIED LOAD VECTOR*/1HO,40X,10HJOINT NO. ,3X,&
5HFX ,5X,5HFY ,9X,5HMZ /)
```

There are some crap characters, sorry. It comes from this thesis from 1972

The book said do it in EXCEL and I was not going to.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

The CDC 6000 series computers used a 48 character set. The quote character was interchangeable with the "not equal" character, ≠ . See PDF page number 28, printed page number I-2-6, in the CDC Fortran Manual .

To add to the confusion, note that the program listing in the thesis appears to have been printed on a console typewriter (IBM Selectric?), which may have substituted an asterisk for the character ≠ . An asterisk in column-1 was also used to indicate comment cards. Thus, there are at least three forms for literal character strings in the program listing.

Regarding the original question, "Is there a faster way to solve this problem in fortran than the long method?", which is ambiguous because the "problem" and the "long method" have not been described, perhaps the HYPOT intrinsic function could be the answer.

`AL(I) = hypot(XM(I), YM(I))`

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

My first experience of FORTRAN was on a CDC Cyber, looking at that pdf manual made me feel quite ill.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Mine was with a CDC 3600, circa 1969.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

I have slowly been fixing the OCR of the code from the Column Program.

This little group is interesting

```
DO 904 I=1, NJL
! READ THE NUMBER OF THE JOINT LOADED, THE LOAD IN THE X-DIRECTION
! THE LOAD IN THE Y-DIRECTION, AND THE CONCENTRATED COUPLE AT THAT JOINT IF ANY.
READ(10,*)JNN(I),FXX(I),FYA(I),FMA(I)
Write(*,55)JNN(I),FXX(I),FYA(I),FMA(I)
Write(12,155)JNN(I),FXX(I),FYA(I),FMA(I)
55 Format(I10,3F10.3)
155 Format(20X,I7,3(3X,F10.3))
904 CONTINUE
DO 1371 I = 1,NJL
JLS= JNN(I)
FYY(JLS) = FYY(JLS) + FYA(I)
FMM(JLS) = FMM(JLS) + FMA(I)
1371 END DO
1336 CONTINUE
```

Notice the use of the FXX and the FYA, why not have an FXA, as the idea is you can load the same joint multiple times.

Bit intrigued. The program will allow up to 1000 load points to the 150 allowed joints.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

FXX is stored in a file on unit 10. It may have use elsewhere than in the loop of 1371.

Jim Dempsey

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

No I checked, the strange idea is he has a set of up to 1000 loaded joints theoretically and then only allocates 100 long vector, setting that aside, the second code sums up for the actual joint number the applied loads. So he could apply 200 units to FYY on Joint 10 and in another line, 300 units to FYY on joint 10 and it would give 500 units, but only the last FXX for the Joint would apply and if it was zero then you have an interesting problem.

He does not publish any results in the thesis, which is a no no, so I cannot compare to anything. But the thesis is worth the read.

I am only getting garbage at the moment, as the goto numbers are in the gutter of the publication and some are missing. I have asked the library to check and told them I will give them the working program so this is not a future problem.

His idea is really good.

His problem in 1972 he did not have instruments that would confirm the results, now we do.

Image if Gauss had had a computer or Turing.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

In the code you provided, FXX is not used (other than to be printed/written). Therefore, I presume, the value was collected, or produced, by a different program that has/had a use for it. The code (to the extent provided), appears to have no computational use for this variable (data point). This is not to say that it is not important to maintain this variable in the dataset.

RE: missing goto numbers

By studying the style of programming, assuming the author was consistent, you may be able to locate receiving points for GOTO (e.g. all may be CONTINUE statements. While this won't explicitly tell you which one, an understanding of what the programmer was coding for will provide you with sufficient information to make an educated guess.

Hopefully you can get your hands on the original publication (or photo copy) as opposed to a mis-scanned .pdf file.

It is unfortunate that the person performing the scanning chose a setup that excluded what were thought to be margins in the document. Hmmm. is that what happened to Fermat's last theorm? (Insignificant scribbling in the margin and of no importance)

Jim Dempsey

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

I wrote to the Librarian at McMaster University and told her of my problem. They have been kind and are going rescan and upload the thesis. I am slowly working through the code with the simple beam problem and it is not hard to work out the kinks. The code is fun to play with - been a while.

You can see the different hands in the code. The 6 by 6 stiffness matrix has been stored in a 3D array of 150 elements, which cuts down the size of the global stiff to 0.3% of the size - saves a lot of zeros.

The code relies on the integer and reals being set to zero when the program starts and there is a lot of assumptions that 0.0 is 0.0 all of the time.

I will upload the code as soon as I am finished, just do not want to waste your time.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Jim:

I ran into an interesting feature on the program, I think that the original developer for the original code back in 1964 to 1967 used a inversion routine that uses the determinant method to solve the inversion. I have not seen it before and the original paper is not found in the library, so I will have to ask.

It is the routine BAND, it was used in a lot of PhD and Masters at McMaster in the period 1965 to 1977. I got as far as band and had to deal with my daughter last night.

It is easier to just replace it, in the beginning with a quick and dirty inverter from Conte and de Boore, but I thought you would like to see the other method.

The A.TXT gives you a summary of the variables up to the point of assembling the global stiffness matrix, my assembly is wrong, I am about to fix it.

I have only just fixed the code the be a minimal run, it is not updated a lot.

Regards

JMN

Again named for a Greek god.

(I have not fixed past band except to get it running with implicit none added. For the moment Real and Integers are enough for the numbers in the program input. )

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

I think you have the GOTO's mixed up. Consider:

```
...
DO 62 J=MP,NM1,M
JP = J - MM
MZC = 0
IF(KK .GE. M) GOTO 1
KK = KK + 1
II = 1
JC = 1
GOTO 2
...
1 Continue
2 Continue
RETURN
```

IOW the only way to get past the DO 62 is if the loop does not execute. Also this makes the code following GOTO 2 through 62 CONTINUE meaningless. So it appears that source you have does not follow the original code (imho).

Jim Dempsey

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Your humble opinion is correct, there are 3 go to numbers that have disappeared into the gutter, I was hoping that it would be obvious, but can I say it is not.

The global stiffness matrix is assembled here

```
DO 521 I = 1,NM
DO 522 JJ = 1,6
IF(MC(I,JJ) .EQ. 0) THEN
Write(12,601)I,JJ,MC(I,JJ)
601 Format(/,20X," Skip joint point Rule A :: ",3(I5))
GOTO 522
end if
DO 523 II =JJ,6
IF(MC(I,II) .EQ. 0) THEN
Write(12,602)I,II,MC(I,II)
602 Format(/,20X," Skip joint point Rule B :: ",3(I5))
GOTO 523
ENDIF
IF(MC(I,JJ)-MC(I,II))524,525,525
525 K=(MC(I,II)-1)*(NB-1) + MC(I,JJ)
A(K)=A(K) + SM(I,JJ,II)
If(K .gt. KMAX) then
KMAX = K
ENDIF
Write(12,603)K
603 Format(/,20X," K value at goto 525 :: ",I5)
GOTO 523
524 K=(MC(I,JJ))*(NB-1)+MC(I,II)
A(K)=A(K) + SM(I,JJ,II)
If(K .gt. KMAX) then
KMAX = K
ENDIF
Write(12,604)K
604 Format(/,20X," K value at goto 523 :: ",I5)
523 END DO
522 END DO
521 END DO
```

I added the 600 lines to provide printout, gives the stiffness matrix as

```
Number of members :: 5
I A(I)
1 0.833E+01
2 0.188E+08
3 0.833E+01
4 -0.978E+06
5 0.469E+08
6 0.833E+01
7 -0.108E+07
8 0.481E+08
9 0.833E+01
10 -0.108E+07
11 0.481E+08
12 0.833E+01
13 -0.108E+07
14 0.481E+08
15 0.293E+08
```

15 elements instead of 225 elements , so either the code is wrong or they have developed a fancy way to only get the non-zeros out of the stiffness matrix. I tend to the first but they could be that clever.

Thanks for the comment, I am assembling the full stiffness matrix and then I will compare it to their 15 numbers and see if I can see sense.

By 1967 surely they had invented inversion routines that did not need the determinant?

Thanks for the help, sometimes you just need someone to tell you, no your not insane, just wrong.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

As you are cleaning up the code, you might as well go all the way and polish it too.

```
DO I = 1,NM
DO JJ = 1,6
IF(MC(I,JJ) .EQ. 0) THEN
Write(12,601)I,JJ,MC(I,JJ)
601 Format(/,20X," Skip joint point Rule A :: ",3(I5))
CYCLE
end if
DO II =JJ,6
IF(MC(I,II) .EQ. 0) THEN
Write(12,602)I,II,MC(I,II)
602 Format(/,20X," Skip joint point Rule B :: ",3(I5))
CYCLE
ENDIF
IF(MC(I,JJ)-MC(I,II) < 0)THEN
K=(MC(I,JJ))*(NB-1)+MC(I,II)
A(K)=A(K) + SM(I,JJ,II)
If(K .gt. KMAX) then
KMAX = K
ENDIF
Write(12,604)K
604 Format(/,20X," K value at goto 523 :: ",I5)
CYCLE
ENDIF
K=(MC(I,II)-1)*(NB-1) + MC(I,JJ)
A(K)=A(K) + SM(I,JJ,II)
If(K .gt. KMAX) then
KMAX = K
ENDIF
Write(12,603)K
603 Format(/,20X," K value at goto 525 :: ",I5)
END DO
END DO
END DO
```

Jim Dempsey

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

How far does polishing go? Once you start it goes on forever! I personally hate all those numbered formats that break up the indentation scheme and make it slower to scan the eye over the code and understand the flow.... I pretty much always have formats as a character variable or parameter statement (my preference).

For scanning the code logic bye eye I don't have IF / ENDIF blocks that only have 1 line in them, 3 code lines rather then 1 seems a waste.....

```
character(*), parameter :: g601 = '(/,20X,A,3(I5))'
DO I = 1, NM
DO JJ = 1, 6
IF( MC(I,JJ) .EQ. 0 ) THEN
Write(12,g601) " Skip joint point Rule A :: ", I, JJ, MC(I,JJ)
CYCLE
end if
DO II = JJ, 6
IF( MC(I,II) .EQ. 0 ) THEN
Write(12,g601) " Skip joint point Rule B :: ", I, II, MC(I,II)
CYCLE
ELSEIF( MC(I,JJ)-MC(I,II) < 0 )THEN
K = (MC(I,JJ))*(NB-1)+MC(I,II)
A(K) = A(K) + SM(I,JJ,II)
If(K .gt. KMAX) KMAX = K
Write(12,g601) " K value at goto 523 :: ", K
CYCLE
ENDIF
K = (MC(I,II)-1)*(NB-1) + MC(I,JJ)
A(K) = A(K) + SM(I,JJ,II)
If(K .gt. KMAX) KMAX = K
Write(12,g601) " K value at goto 525 :: ", K
END DO
END DO
END DO
```

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Jim:

I am struggling against an internal deadline, so the code clean up can wait a while. I got to the BMPCAL routine, it calculates a new axial load and bending moment for the column, took me a few hours with the code and a pencil and the TV mumbling in the background to work out where the go to statements went to. 100 was in the wrong place as was 30. Now it gives answers I understand.

Now to work out MPHI.

Have you seen this before?

If so, have you tried it?

The real problem is this Canadian team was on the right track, but they did not have the instruments to do the work on real bridges, not till the accurate accelerometer came along. They did the work with compression machines with dial gauges and hand control. They got a long way and then ran into the wall. Now everyone wants to use strain gauges and SSI, it is hard to explain to them that Fortran is the solution not MAPLE etc...

How do I clean up the goto 4, this takes the full assembled stiffness matrix and then eliminates the restrained elements, allows me to check full matrix, temp holds the rows and columns to be eliminated.

This program will never handle more than 10 members as it is a checking program. KAMAX is the reduced size of matrix and LM1 the number to eliminate. LM = LM1+KAMAX

```
I1 = 1
J1 = 1
Do 3 I = 1,KAMAX
do 4 J = 1,KAMAX
do 7 K1= 1,LM1
If(I .eq. TEMP(K1) .or. J .eq. TEMP(K1)) then
!Write(*,*)I,J,"SKIPPED"
GOTO 4
endif
7 END DO
! Write(*,*)I,J,I1,J1
SKR(I1,J1) = Sk(I,J)
IF(J1 .eq. LM) then
J1 = 1
I1 = I1+1
else
J1 = J1 + 1
end if
4 end do
3 end do
```

Andrew: Neat idea.

```
IF(W2 .EQ. 0.0) W2 = 1.0
IF(W3 .EQ. 0.0) W3 = 1.0
WA = ABS(W2)
WB = ABS(W3)
TEMP = W2*(WA+WY-ABS(WA-WY) )/(2.0*WY*WA)
STEEL2= FSY*TEMP
STEEL3= FSY*W3*(WB+WY-ABS(WB-WY) )/(2.0*WY*WB)
IF(W2 .EQ. 0.0) STEEL2 = 0.0
IF(W3 .EQ. 0.0) STEEL3 = 0.0
```

W2 and W3 are never zero, I had to chuckle.

@mecej4 - thanks again, I used the LINSOLVER you created to take out the strange inverter. No idea how their code worked, much quicker and easier and simple to check results.

JMN

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

I've had issues with Edit and Continue on IVF before (long time ago) and have just turned it off. The Build process is quite fast (assuming you don't change something that every module uses).

>>How do I clean up the goto 4

You name the loop, then issue CYCLE TheNameYouUsed

But this is just a style change, not addressing the issue mecej4 raises.

```
I1 = 1
J1 = 1
loop_I: Do I = 1,KAMAX
loop_J: do J = 1,KAMAX
do K1= 1,LM1
If(I .eq. TEMP(K1) .or. J .eq. TEMP(K1)) then
!Write(*,*)I,J,"SKIPPED"
CYCLE loop_J
endif
END DO
! Write(*,*)I,J,I1,J1
SKR(I1,J1) = Sk(I,J)
IF(J1 .eq. LM) then
J1 = 1
I1 = I1+1
else
J1 = J1 + 1
end if
end do
end do
```

Jim

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Use of FINDLOC can clean it up too.

```
I1 = 1
J1 = 1
Do I = 1,KAMAX
if(FINDLOC(TEMP(1:LM1), VALUE=I) /= 0) CYCLE
do J = 1,KAMAX
if(FINDLOC(TEMP(1:LM1), VALUE=J) /= 0) CYCLE
! Write(*,*)I,J,I1,J1
SKR(I1,J1) = Sk(I,J)
IF(J1 .eq. LM) then
J1 = 1
I1 = I1+1
else
J1 = J1 + 1
end if
end do
end do
```

Jim Dempsey

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page