Get rid of goto jump

[JohnNichols]

IF(ISW .eq. 1 )then
                ID=0
                M=0
                L=NDF

                !-------------------------------------------------------------------------------------------------------
                !
                !  
                !
                !------------------------------------------------------------------------------------------------------

903             JJ=J1-1
                IF(JJ .gt. 0) then
905                 DO 910 J=1,JJ
                        DO 910 K=1,6
910                 ID = ID + JTYPE(J,K)
                endif
912             DO 965 IC=1,6
                    IF(JTYPE(J1,IC) .lt. 1)then
915                     DO 916 N=1,NLS
916                     ASAT(M + IC,N) = 0.0
                    else
920                     DO 925 N=1,NLS
925                     ASAT(M+IC,N)=GK(ID+1,L+N)
                        ID=ID+1
                    endif
965             end do
                IF(M .lt. 1) then
980                 M=6
                    J3=J1
                    J1=J2
                    ID=0
                    GO TO 903
                endif

                !-------------------------------------------------------------------------------------------------------
                !
                !  
                !
                !---------------------------------------------------------------------------------------------------

990             DO 1100 II=1,6
                    DO 1100 J=1,NLS
                        SR(II,J)=0.
                        DO 1100 K=1,12
1100            SR(II,J)=SR(II,J)+SAT(II,K)*ASAT(K,J)
                !Write(*,1010)I
                if(i .eq. 1) then
                    Write(*,1011)LSN(I)
                    Write(sw,1011)LSN(I)
                endif
                DO 1020 J=1,NLS
                    write(*,1030)i,(SR(II,J),II=1,6),(-(SR(2,j)+SR(3,j))/Length(I)),(-(SR(4,j)+SR(4,j))/Length(I))
1020            write(sw,1030)i,(SR(II,J),II=1,6),(-(SR(2,j)+SR(3,j))/Length(I)),(-(SR(4,j)+SR(4,j))/Length(I))
                !pause
1030            FORMAT(I4, 5x,1X,5(F11.4,2x),F10.4,1x, F10.4, 1x, F10.4)
                N=4
                J=12
            endif

I have been getting rid of all the arithmetic goto's inside this "old" Fortran -- from the 1973 book from Harrison --  I am blowed if I can see an elegant method to get rid of the goto 903. I have had several shots at it and I end up always breaking the program. 

Any ideas would be appreciated - I was thinking a do while on M?

[andrew_4619]

If it works leave it alone is my advice! If you need to tinker why not get rid of most of the Do loop by using array opporations it will make it easier to read/understand ......

[mecej4]

I second andrew_4619's advice. In fact, ditch the old spaghetti and rewrite the whole code, if you want modern style and understand the algorithm completely. If you tinker with old code with lots of GOTO and IF statements, it is very easy to introduce new bugs, unless you use a proven code polisher.

One idea: lines 12-34 could be replaced with (please double check; not listening to my own advice here!):

do
   do j=1,j1-1
      id=id+sum(jtype(j,1:6))
   end do
   do ic=1,6
      if (jtype(j1,ic) < 1) then
         asat(m+ic,1:nls) = 0
      else
         asat(m+ic,1:nls)=gk(id+1,l+(1:nls))
         id=id+1
      endif
   end do
   if (m >= 1) exit
   m=6
   j3=j1
   j1=j2
   id=0
end do 

[FortranFan · Accepted solution]

John Nichols wrote:

.. I have been getting rid of all the arithmetic goto's inside this "old" Fortran -- from the 1973 book from Harrison --  I am blowed if I can see an elegant method to get rid of the goto 903. I have had several shots at it and I end up always breaking the program. 

If your arrays are not too large for you to consider vectorization/parallelization options, you can consider starting from something as shown below: naturally UNTESTED and with the two "magic numbers" of 6 and 12 replaced by named constants of MN_6 and MN_12 respectively.

   if ( ISW == 1 ) then

      do K = 1, 2
         
         ID = 0
         M = (K-1)*MN_6

         if ( J1 > 1 ) then
            ID = sum( JTYPE( 1:J1-1, 1:MN_6) )
         end if
         
         do IC = 1, MN_6
            if ( JTYPE(J1,IC) < 1 ) then
               ASAT(M+IC,1:NLS) = 0.0
            else
               ASAT(M+IC,1:NLS) = GK(ID+1,NDF+1:NDF+NLS)
               ID = ID + 1
            end if
         end do
         
         if ( K == 2) then
            J3 = J1
            J1 = J2
         end if
   
      end do      

      SR(1:MN_6,1:NLS) = matmul( SAT(1:MN_6,1:MN_12), ASAT(1:MN_12,1:NLS) )
      
      if ( i == 1 ) then
         write(*,1011)  LSN(I)
         write(sw,1011) LSN(I)
      end if
      
      do J = 1, NLS
         write(*,1030) i,(SR(II,J),II=1,6),(-(SR(2,j)+SR(3,j))/Length(I)),(-(SR(4,j)+SR(4,j))/Length(I))
         write(sw,1030)i,(SR(II,J),II=1,6),(-(SR(2,j)+SR(3,j))/Length(I)),(-(SR(4,j)+SR(4,j))/Length(I))
1030     format(I4, 5x,1X,5(F11.4,2x),F10.4,1x, F10.4, 1x, F10.4)
      end do
      
      !pause
      N=4
      J=12

   end if

P.S.> Code edited: added matmul in place of dot_product in the assignment statement for SR.

[JohnNichols]

Dear All:

I cannot help tinkering - I am an engineer.  You are not on this forum unless you enjoy some form of pain - daily.

I understand the algorithm very well, studied it for 3 year at Uni and I can do it using an plain inverter, but this program came with a nice set of pictures in the book and some samples.  It will be rewritten in entirely new form by experts, I merely hack a crude path in the wilderness to see what is possible.  An eigensolver has been fitted to the rear end and it is providing some light relief.  MKL  Feast is excellent --

Thanks for the help I will now try them.

Regards

John

[JohnNichols]

asat(m+ic,1:nls)=gk(id+1,l+(1:nls))

 The complier did not like this or the dot product - thanks for trying.

[andrew_4619]

asat(m+ic,1:nls)=gk(id+1,L+1:L+nls)

looks righter...

[TimP]

You could write that loop beginning at label 903 as

do m=0,6,6

....

  ID=0

  end do

I also have a personal preference for sanitizing code this way.  Why make it unnecessarily obscure?

The if(jj > 0) was made redundant by F77, although F66 didn't define that case without the explicit comparison.

MATMUL is the rank 2 equivalent of dot_product.

[JohnNichols]

Thanks

Matmul worked and I got the standard answers as before so the replacement worked - again thanks. Still have to do a full check - but it appears to work.

I would post the code - but a very expensive lawyer explained IP to me and made me aware of the dangers of publishing code that has a value

Interestingly the Eigen solver does not give exact multiples for the harmonics -- off by a little bit.  The MKL solver is way faster than the eigensolvers in the structures program and it gives the vectors which really help.

John