Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.
28435 Discussions

## Optimal Code

Valued Contributor III
2,419 Views
``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?

24 Replies
Honored Contributor III
2,130 Views

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

Jim Dempsey

Valued Contributor III
2,116 Views

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

Honored Contributor III
2,106 Views

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

Spaces matter.

Jim Dempsey

Valued Contributor III
2,097 Views

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

Honored Contributor III
2,085 Views

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))``
Honored Contributor II
2,078 Views

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

Honored Contributor III
2,070 Views

Mine was with a CDC 3600, circa 1969.

Valued Contributor III
1,998 Views

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

This little group is interesting

``````DO  904 I=1, NJL

!   THE LOAD IN THE Y-DIRECTION, AND THE CONCENTRATED COUPLE AT THAT JOINT IF ANY.

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.

Honored Contributor III
1,997 Views

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

Jim Dempsey

Valued Contributor III
1,987 Views

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.

Honored Contributor III
1,970 Views

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

Valued Contributor III
1,961 Views

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.

Valued Contributor III
1,930 Views

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

Honored Contributor III
1,919 Views

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

Valued Contributor III
1,915 Views

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.

Honored Contributor III
1,863 Views

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

Honored Contributor II
1,845 Views

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

Valued Contributor III
1,779 Views

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

Honored Contributor III
1,738 Views

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

Honored Contributor III
1,738 Views

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