Numerical Recipes in Fortran Gaussian elimination

[JohnNichols]

NR does not return the correct inverse for the test case on Gaussian elimination with current IFX. 

I do not have time to track down the error, it may be something I did, but this is an exact copy of the original with the same matrix. 

It gives the original matrix back.  

Plus I have been denied access to the site all weekend. 

[JohnNichols]

Jim:

I was playing with some old code from Harrison from 1973, he provides a second order elastic analysis, I want to add it to ULARC, spent four days on the code, it has taken me all that time to realize he had a I in a do loop instead of a 1, why would we use I as the common counting variable.   do J = I,N is so easy to miss in a lot of code.

The interesting feature and I will shoot the code along as soon as I clean it up, is the first answer is a 3 by 3 array for the deflections but element 2,2 does not match the IFX answer, he has -2.9 and I have -5.7,  all the others are about the same, some to 5 places.  

I think he had a typo, I found a couple, I wonder how many typos exist in current code?  

Using IFX instead of IFORT halves the run time.  The other interesting feature, the code lacked an invertor for the stiffness matrix.  He assumes you are smart enough to know that and steal one from the earlier code in the book.  

He uses a different method to establish the stiffness matrix, he pushes the elements into two vectors and then builds the stiffness matrix from the vectors, it saves some arithmetic steps.  Not seen that before.  

I did not want to rip all of mecej4's code out of the 3D solver for MKL inversion, so I was going to use the NR, but it would not invert with IFX, so I use the old Conte and de Boore.  

Anyway I hope all is well at your place.  

John

[JohnNichols]

Book with code.  Flood limit is gone.

[andrew_4619]

"do J = I,N is so easy to miss in a lot of code." IMPLICIT NONE not used maybe? A good tool for exposing typo's

[JohnNichols]

This was not a typo, the code is looking to review all of the off diagonal elements, by my missing that I had a few hours of pain as I tried to understand the code and then I looked at the original again, as it just did not make sense.  It works now.  Implicit none will not pick up the difference between a 1 and I in a do loop, when the original author uses I for almost every do loop. 

Steve:

My apologies, I do not want to waste your time, it occurred several times when I was compiling the code over a short time, then it stopped, I was just intrigued, I was wondering if it was a result of deleting the format, and then hitting debug, maybe the compiler missed the delete or it did not rebuild.  It is not something to worry about, just oddish. 

NR, my mistake, I should have spotted that.  

Thanks for all the comments. 

John

Attached is the code from Harrison for a second order analysis, this is the last bit to sort out - remove the arithmetic ifs, and work out the algorithm, although it runs nicely.  

            ZZ = 0.00000005
            flag = 0
            TESTA = ZERO
            Do 700 I = 1,L

                  !!TESTA = TestA + DUM
                  !if(I .eq. L) then
                  ! write(*,*)I,TestA
                  ! endif
                  !   OLDDEFX(I) = DEFX(I)
                  Test = ABS(DEFX(I)*ZZ)
                  If(ABS(DEFX(I) - DELTA))700,690,690
690               IF(ABS(DEFX(I) - ASAT(I,KJ) - TEST))700,700,770
                  ! if(test < delta) then
                  !  goto 950
                  ! end if

700         end do
            goto 950

[avinashs]

@JohnNichols The error is due to the fact that you are using the Fortran 77 driver program for the F90 version of Gaussian elimination. The code works fine if you use the F90 driver. In the F77 driver, array dimensions are explicitly passed whereas the F90 version has assumed shaped arrays so all arrays must be allocated with exact dimensions (n,m). Within GAUSSJ, the F90 version determines (n,m) using the size operator on A and does not pass (n,m) explicitly. I had checked the NR F90 codes in year 2000 and found all to be correct except SVD, which was later corrected.

[andrew_4619]

well the ground base of investigation was to add  warn:all i.e. ask the compiler to look for problems rather than ignore them. Et Voila!