Help with some old Fortran code

[BridgeGuy]

I need help understanding how this code flows. Back when I was in school, we had an assignment to create a 3-d truss structural analysis program. It became a colaborative effort as most of us had little programming experience. We came across this subroutine which I believe performs a Gauss elimination and back substitution. I am currently trying to teach myself how to program in Python and I'm attempting to translate the old Fortran code into Python 3.x. Here is the subroutine.

 SUBROUTINE BNDSOL(BGK,Q,NDIS,MB)
      DIMENSION BGK(100,30),Q(100),F(30)
      N=0
  500 N=N+1
      Q(N)=Q(N)/BGK(N,1)
      IF(N-NDIS)550,700,550
  550 DO 600 K=2,MB
      F(K)=BGK(N,K)
  600 BGK(N,K)=BGK(N,K)/BGK(N,1)
      DO 660 L=2,MB
      I=N+L-1
      IF(NDIS-I)660,640,640
  640 J=0
      DO 650 K=L,MB
      J=J+1
  650 BGK(I,J)=BGK(I,J)-F(L)*BGK(N,K)
      Q(I)=Q(I)-F(L)*Q(N)
  660 CONTINUE
      GO TO 500
  700 N=N-1
      IF(N)750,900,750
  750 DO 800 K=2,MB
      L=N+K-1
      IF(NDIS-L)800,770,770
  770 Q(N)=Q(N)-BGK(N,K)*Q(L)
  800 CONTINUE
      GO TO 700
  900 RETURN
      END

BGK is the structure stiffness matrix, Q is the applies loads, NDIS is the number of degrees of freedom and MB is the bandwidth of BGK. As I understand this, these four variables are passed to the subroutine and BGK and Q is passed back. I believe the code up to line 700 is the elimination portion of the routine and from 700 on is the back substitution. Please correct me if I'm wrong. How I understand the code is that it begins, increments N to 1 and divides the applied load in row 1 by the pivot value. It then checks to see if N-NDIS is non-zero. if it isn't it runs a do loop to 600 creating the vector F from BGK and dividing the rest of row 1 by the pivot value. The next do loop to 660 begins by performing another check. If the check is >= to 0 then it runs another do loop to 650 or breaks out of the do loop if the number is negative. Once the nested do loop to 650 finishes and the do loop to 660 completes, the routine returns to line 500 and increments N and the loops restart. The go to 500 line seems to function like a for statement. Lines 700 to 900 operate similar to the first part of the routine except it increments N by -1. Have I correctly interpreted the subroutine? I appreciate any help I can get.

 

[xwb]

First translate from Fortran II to F77, incorporating your comments so that people who can't read FII can understand it. Might be easier for you to translate to Python in this form.

      SUBROUTINE BNDSOL(BGK,Q,NDIS,MB)
      DIMENSION BGK(100,30) ! stiffness matrix
      DIMENSION Q(100)      ! loads
      INTEGER NDIS          ! degrees of freedom
      INTEGER MB            ! bandwidth of BGK
      DIMENSION F(30)
      ! Possibly elimination
      N=0
      DO
         N=N+1
         Q(N)=Q(N)/BGK(N,1) ! not sure about this
         IF(N .EQ. NDIS) EXIT
         ! Divide by pivot value
         DO  K=2,MB
            F(K)=BGK(N,K)
            BGK(N,K)=BGK(N,K)/BGK(N,1)
         END DO
         DO  L=2,MB
            I=N+L-1
            IF(NDIS .GE. 0) THEN
               J=0
               DO K=L,MB
                  J=J+1
                  BGK(I,J)=BGK(I,J)-F(L)*BGK(N,K)
               END DO
               Q(I)=Q(I)-F(L)*Q(N)
            ENDIF
         END DO
      END DO
      ! Possibly back substitution
      DO N = N-1, 0, -1
         DO K=2,MB
            L=N+K-1
            IF(NDIS .GE. L) THEN
               Q(N)=Q(N)-BGK(N,K)*Q(L)
            END IF
         END DO
      END DO

      RETURN
      END

 

[xwb]

The do loops should be

      DO N = 1, NDIS-1, 1
         Q(N)=Q(N)/BGK(N,1) ! not sure about this
         ...
      END DO
      Q(NDIS) = Q(NDIS) / BGK(NDIS, 1)
      DO N = NDIS - 1, 1, -1
         ...
      END DO

 

[salgerman]

Hhhhmmmm...I have not looked at the code above; but, with Python being so rich in libraries and, for sure, with numpy and scipy, why translate a matrix solver? Python has a few of its own.

O.k., I couldn't resist, I peeked at the first line...why is BGK not square? Since stiffness matrices are symmetric (I just looked it up), this non-square storage gives me a hint that they are only storing the upper triangular portion of the stiffness matrix and that they have left-justified every row so that the diagonal values are all on column 1.

SO!...Before you re-use this piece of code in Fortran OR Python, you'd better know exactly what it does and what arguments is expects.

gsal

 

[BridgeGuy]

Thanks for the replies and suggestions. For the example I'm using to hand check the main body of code, the stiffness matrix is square and with only the upper triangle. The elements of the matrix are not left justified though.
The methodology of the main program combines the individual element stiffness matrices into the structure stiffness matrix using the code number technique. If I remember correctly, this technique ensures a square upper triangle matrix. I'm not sure why the dimensions in the subroutine are not square. I didn't work on this part of the program as another person on our team did the programming for the subroutine.
I agree with gsal about using one of python's libraries from numpy or scipy but I wanted to make sure I understood what was happening in this subroutine first. If it would help clarify how the subroutine is used, I can post the full program code. I appreciate everyone's time and suggestions.

 

[salgerman]

BridgeGuy:

You say that the elements of the matrix are not left justified...I say they are, but maybe it was not clear what I meant by 'left-justified' ...if you open a text file with a fully populated square stiffness matrix in it and from every row you use the delete key to delete every value from every column up to, but no including, the diagonal value...that's what I meant by left-justified...and so, the diagonal values are always in column 1 of the new storage scheme.

You can see in the algorithm of the routine that they divide by BGK(N,1) in a couple of places...why do you think they do that? what's special about the value in (N,1)?

Then, if you go back to the text file and you delete the largest rectangular area of zeroes...now you have the needed storage and the number of columns in this new "rectangular" matrix is the bandwidth.

The reason why it is important to know what they are doing and HOW the are handling BGK is because if you pass a full matrix or even if you pass a square matrix that is only populated on the upper triangular...that subroutine is not going to work correctly...for it to work correctly you need to pass the stiffness matrix stored in the expected manner.

Exercise: 9x9 stiffness matrix

Full matrix

q w e 0 0 0 0 0 0 
w s d f 0 0 0 0 0 
e d u j 0 0 0 0 0
0 f j c v b 0 0 0
0 0 0 v h w e 0 0
0 0 0 b w u h y 0
0 0 0 0 e h r e r 
0 0 0 0 0 y e w d 
0 0 0 0 0 0 r d j

Remove columns from row up to but no including diagonal

q w e 0 0 0 0 0 0 
s d f 0 0 0 0 0 
u j 0 0 0 0 0
c v b 0 0 0
h w e 0 0
u h y 0
r e r 
w d 
j
[code]

Remove zeros on the right of all rows starting from the one on the right of the right-most non-zero from all rows
[code]
q w e 
s d f 
u j 0
c v b
h w e
u h y
r e r 
w d 
j

Oh...just found this, check it out and page 8/15...if I had looked for it earlier, I wouldn't have had to type all this ;-)