Optimal Code

JohnNichols · ‎03-21-2023

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?

jimdempseyatthecove · ‎03-21-2023

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

JohnNichols · ‎03-21-2023

*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 " "

jimdempseyatthecove · ‎03-21-2023

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

Spaces matter.

Jim Dempsey

JohnNichols · ‎03-21-2023

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.

mecej4 · ‎03-21-2023

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

andrew_4619 · ‎03-21-2023

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

mecej4 · ‎03-21-2023

Mine was with a CDC 3600, circa 1969.

JohnNichols · ‎03-23-2023

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.

jimdempseyatthecove · ‎03-23-2023

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

Jim Dempsey

JohnNichols · ‎03-23-2023

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.

jimdempseyatthecove · ‎03-23-2023

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

JohnNichols · ‎03-23-2023

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.

JohnNichols · ‎03-24-2023

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

jimdempseyatthecove · ‎03-24-2023

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

JohnNichols · ‎03-24-2023

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.

jimdempseyatthecove · ‎03-25-2023

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

andrew_4619 · ‎03-25-2023

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

JohnNichols · ‎03-26-2023

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

jimdempseyatthecove · ‎03-27-2023

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

jimdempseyatthecove · ‎03-27-2023

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