Part 11. A Complete Fortran Program

(c) 2017 by Barton Paul Levenson

A Gauss-Jordan elimination program

This is a full-scale Fortran program that actually does something useful. It performs Gauss-Jordan elimination on a matrix in order to solve a system of linear equations. If you don't know what that means, see Appendix 4 of the tutorial on statistics.

The basic code

Here is a module to hold the global variables:

module matrixGlobals
    implicit none
    integer, parameter :: d = selected_real_kind(p = 15)
    integer, parameter :: maxRows = 20

    real(d) :: aug(maxRows, maxRows + 1)	! Augmented matrix
    integer :: Nrows, Ncols			! Number of rows, columns.
end module matrixGlobals

The variable aug is the matrix itself. It has one more column than it has rows, the column for the answers.

We will use a top-down approach to write this program. That means breaking it into a separate routine for each clearly defined function. Major variables have to be shared across routine boundaries; thus the module.

Here's the main program:

! main program
program matrix
    call setup				! Load the matrix.
    call showMat			! Check that we got it right.
    call Gauss				! Do the matrix math.
    call showMat			! Check the results.
end program matrix

The logic is straightforward. We initialize the matrix. We check that we did so properly. We perform Gauss-Jordan elimination. Then we look at the transformed matrix. If everything has gone right, the answers will be in the rightmost column.

Before we write any of these functions, we need to know what data we will be working on. Let's start small, with a 3 by 4 matrix. We will have three variables, X, Y, and Z, which we want to find the values of. Let's put some totally arbitrary constants in front of them and leave the column of constants empty:

Now, let's see what answers we want, since we have to know in advance to check if the program works properly. Say x should be 1, Y 2, and Z 3. Plugging those in, we find:

So at the beginning, our matrix should look like this:

If our program works correctly, the output should look like this:

Here's the code for the setup routine:

! setup initializes the matrix.
subroutine setup
    use matrixGlobals
    implicit none

    aug(1, 1) = 1d0; aug(1, 2) = 2d0; aug(1, 3) = 3d0; aug(1, 4) = 14d0
    aug(2, 1) = 4d0; aug(2, 2) = 9d0; aug(2, 3) = 8d0; aug(2, 4) = 46d0
    aug(3, 1) = 5d0; aug(3, 2) = 6d0; aug(3, 3) = 1d0; aug(3, 4) = 20d0

    Nrows = 3
    Ncols = Nrows + 1
end subroutine setup

It ain't pretty, but it works.

To check whether we put in the proper values in the proper order, we need a routine to print out what's in the matrix:

! showMat displays what's in the matrix.
subroutine showMat
    use matrixGlobals
    implicit none
    integer  :: i, j					! Loop counters.

    write (*, *)					! Print a blank line.

    do i = 1, Nrows
        do j = 1, Ncols
        	write (*, '(f7.1, $)') aug(i, j)
        end do

        write (*, *)					! Go to next line.
    end do

    write (*, *)					! Print one last blank line.
end subroutine showMat

Note the "nested loop" where one do loop indexed by i encloses another, indexed by j. Nested loops are very common when programming with two-dimensional arrays. Note how j counts through its entire range of values before i gets to its next value.

At this point, you might try commenting out the call to Gauss, and the second call to ShowMat, in the main program, and running what we have so far. If all goes well, you should see this:

We need one more major function to make the whole thing run properly--the Gauss function.

Gauss-Jordan elimination

Let's go over how we would solve this problem by hand. The matrix is:

We work our way through it one column at a time, but only with as many columns as there are rows--we concentrate on the "square matrix" at the left.

For the first column, the element on the main diagonal is "1," so there's no need to scale this row. But we need zeroes in the other column elements, so we have to scale the 2nd and 3rd rows. That will put 1s where we now have 4 and 5.

Once we've done that, we subtract each row from the top row to get zeroes in the required places. This operation is called "sweeping" the column, as if we were sweeping away the non-zero values.

Here it is a bit at a time. First, we scale the 2nd row by multiplying through by 1/4:

For now, we won't bother reducing 46/4 to 23/2.

First row minus second row replaces the second row. This is "row subtraction:"

Now we scale the third row by 1/5:

And we do row subtraction, 1st row minus 3rd:

Now we start on the second column. The diagonal element is now -1/4 rather than 1, so we must scale that row: multiply it by -4:

Then we scale rows 1 and 3 and do row subtraction from row 2, getting:

We've swept the column properly, but as a result, we have a problem: The diagonal element in the first column is no longer 1. In fact, if the whole system is swept, we get:

So the final step must be to rescale each row to get the "identity matrix" on the left--all 1s on the main diagonal, 0s elsewhere:

Gauss-Jordan elimination code

We can write the Gauss subroutine as follows. Notice we have deferred the scaling and sweeping operations to subordinate subroutines which Gauss calls.

! Gauss solves a system of linear equations by Gauss-Jordan
! elimination.
subroutine Gauss
    use matrixGlobals
    implicit none

    integer  :: col				! Column loop counter.
    real(d)  :: scaleFactor		! Scale factor.

    do col = 1, Nrows
        call sweep(col)
    end do

    do col = 1, Nrows
        scaleFactor = 1d0 / aug(col, col)
        call scale(col, scaleFactor)
    end do
end subroutine Gauss

First we sweep each column. Then we rescale the columns.

In that last loop, the scale factor is always the reciprocal of the element on the main diagonal. The elements on the main diagonal are A1,1; A2,2; A3,3; and so on if we have a larger matrix. Thus we only need to know the column number. By definition, on the main diagonal, the row and column numbers are always equal.

Subordinate functions

scale

! scale multiplies a row by a constant.
subroutine scale(row, factor)
    use matrixGlobals
    implicit none

    integer, intent(in) :: row			! Row number.
    real(d), intent(in) :: factor		! Scale factor.

    aug(row, 1:Ncols) = aug(row, 1:Ncols) * factor
end subroutine scale

In most languages we would need to loop across columns to do this. In modern Fortran, we can use the array-handling capabilities built into the language, and specify a range of indices (1:Ncols) in the array in question.

sweep

! sweep puts zeroes in all elements of a column but the pivot.
subroutine sweep(pivot)
    use matrixGlobals
    implicit none

    integer, intent(in) :: pivot		! Base row number.
    real(d) :: factor			! Scale factor.
    integer :: row, col			! Row, column loop counters.

    factor = 1d0 / aug(pivot, pivot)	! Set the scale factor.
    call scale(pivot, factor)		! Scale the base row.

    do row = 1, Nrows
        if (row .ne. pivot) then
            factor = 1d0 / aug(row, pivot)
            call scale(row, factor)

            do col = 1, Ncols
                aug(row, col) = aug(pivot, col) - aug(row, col)
            end do
        end if
    end do
end subroutine sweep

This looks intimidating, but if you look it over slowly, a line at a time, you can see what it's doing. The first two lines merely scale the row we're working from--the "pivot row" or "base row," indexed by the variable pivot. That makes the element on the main diagonal 1.

Next comes the big nested structure. Outermost is a do loop, going from the first row to the last. Inside that is an if structure. The inside of the loop executes only if we are not in the pivot row. We have a 1 there where we want it.

Inside the if structure is the heart of the function:

            factor = 1d0 / aug(row, pivot)
            call scale(row, factor)

            do col = 1, Ncols
                aug(row, col) = aug(pivot, col) - aug(row, col)
            end do

First the target row is scaled to 1, so row subtraction will yield a 0 in the swept element. The loop simply goes through the target row, one column at a time, and subtracts the element from the pivot row, putting the answer where the target element was. When all that is done for each row, we're done.

Now we have a whole program: data module, main routine, plus the subroutines scale, setup, Gauss, showMat, and sweep.

Write it and run it, and we get:

The setup routine may be awkward if you want to do a 20 by 20 matrix. Can you modify the program to read the matrix elements from a file?

And what if the matrix is singular? Suppose you have a 0 on the main diagonal? Can you bulletproof the program against such cases, or at least trigger a warning message to the user?