! This is a Fortran 90 test program; the data are the same as used in
! the previous example, converted to Fortran syntax.  The use of
! TRANSPOSE is because arrays are stored the other way round to
! Mathematica.

! It uses only whole-row operations, whereas most descriptions use
! only the sections of the rows that have not already been zeroed.
! This is only because it is a little simpler - the other form would
! work just as well.

! Essentially, I got Gaussian Elimination working in this, converted it
! to the equivalent Mathematica, and THEN added the extras to make it
! handle symbolic values.  People who don't know Fortran would find it
! easier to develop the simple, floating-point Mathematica straight off.
! Do whatever is easiest.

! Note that I have left some debugging PRINT statements in.

PROGRAM Main
    IMPLICIT NONE
    INTEGER, PARAMETER :: DP = KIND(0.0D0), n = 5
    REAL(KIND=DP) :: A(N,N) = RESHAPE( (/ &
           0.483836, 0.266372, 0.579677, 0.234781, 0.0530176,  &
           0.759127, 0.575231, 0.547177, 0.77256, 0.0687799,  &
           0.649651 , 0.942712 , 0.0837541 , 0.186192 , 0.235544 , &
           0.236206, 0.292983, 0.349548, 0.290192, 0.475288,  &
           0.522635, 0.863739, 0.689245, 0.755017, 0.997574 &
       /),(/5,5/) )
    REAL(KIND=DP) :: C(N) = (/ &
           0.203949, 0.887397, 0.4687, 0.624823, 0.200951 &
       /)
    REAL :: x
    INTEGER :: i, j
    A = TRANSPOSE(A)
    PRINT '(5(5F10.6/))', TRANSPOSE(A)
    PRINT '(5F10.6/)', C
    DO j = 1,N-1
        DO i = j+1,N 
            x = A(i,j)/A(j,j)
            A(i,:) = A(i,:)-x*A(j,:)
            C(i) = B(i)-x*B(j)
        END DO
    END DO
    PRINT '(5(5F10.6/))', TRANSPOSE(A)
    PRINT '(5F10.6/)', C
    DO j = N,2,-1
        DO i = 1,j-1
            x = A(i,j)/A(j,j)
            A(i,:) = A(i,:)-x*A(j,:)
            C(i) = B(i)-x*B(j)
        END DO
    END DO
    PRINT '(5(5F10.6/))', TRANSPOSE(A)
    DO i = 1,N
        C(i) = B(i)/A(i,i)
    END DO
    PRINT '(5F10.6/)', C
END PROGRAM Main