×
INTELLIGENT WORK FORUMS
FOR COMPUTER PROFESSIONALS

Are you a
Computer / IT professional?
Join Tek-Tips Forums!
• Talk With Other Members
• Be Notified Of Responses
• Keyword Search
Favorite Forums
• Automated Signatures
• Best Of All, It's Free!

*Tek-Tips's functionality depends on members receiving e-mail. By joining you are opting in to receive e-mail.

#### Posting Guidelines

Promoting, selling, recruiting, coursework and thesis posting is forbidden.

# Help with some old Fortran code

## Help with some old Fortran code

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

#### CODE --> Fortran

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.

### RE: Help with some old Fortran code

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.

#### CODE

SUBROUTINE BNDSOL(BGK,Q,NDIS,MB)
DIMENSION BGK(100,30) ! stiffness matrix
INTEGER NDIS          ! degrees of freedom
INTEGER MB            ! bandwidth of BGK
DIMENSION F(30)
! Possibly elimination
N=0
DO
N=N+1
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 

### RE: Help with some old Fortran code

The do loops should be

#### CODE

DO N = 1, NDIS-1, 1
...
END DO
Q(NDIS) = Q(NDIS) / BGK(NDIS, 1)
DO N = NDIS - 1, 1, -1
...
END DO 

### RE: Help with some old Fortran code

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

### RE: Help with some old Fortran code

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

### RE: Help with some old Fortran code

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

#### CODE

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

#### CODE

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

#### Red Flag This Post

Please let us know here why this post is inappropriate. Reasons such as off-topic, duplicates, flames, illegal, vulgar, or students posting their homework.

#### Red Flag Submitted

Thank you for helping keep Tek-Tips Forums free from inappropriate posts.
The Tek-Tips staff will check this out and take appropriate action.

Close Box

# Join Tek-Tips® Today!

Join your peers on the Internet's largest technical computer professional community.
It's easy to join and it's free.

Here's Why Members Love Tek-Tips Forums:

• Talk To Other Members
• Notification Of Responses To Questions
• Favorite Forums One Click Access
• Keyword Search Of All Posts, And More...

Register now while it's still free!