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

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.

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

## CODE

## RE: Help with some old Fortran code

## CODE

## RE: Help with some old Fortran code

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

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

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

Remove columns from row up to but no including diagonal

## CODE

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