decomp is a simple implementation in Fortran of the LU-decomposition from the book Computer Methods for Mathematical Computations. It is just 169 lines and theoretically one can use the code to learn the algorithm. It also allows you to get the flavor of Fortran, the code, I guess, is written in Fortran IV. Fortran is a simple language and it should be possible to learn it, just by reading the code.

What is important to remember is that a matrix could be close to singular (ill-conditioned) and in this case the numerical solution is very sensitive to rounding errors. This implies that it is good to estimate the condition number of the matrix

http://en.wikipedia.org/wiki/Condition_number

to judge the results.

solve is the back substitution subroutine (only 57 lines of code) and it should be used after decomp provided that the condition number estimated by decomp is not too high.

I have used the main program from the book Numerical Methods and Software and just have modified it accordingly to call decomp and solve. The code here is already in Fortran 77. The reason to take the main program from another book is that I am just lazy and did not want to write it by myself from the scratch (I will do it in the next section in C++). I should say that the software to solve a linear system in  Numerical Methods and Software is better but more complex.

The main disadvantage of decomp and solve is that they do not use BLAS and their performance will not be that good. Yet, it is a good reference point to check what improvements an optimized BLAS will bring.

Please look at the code of decomp, solve and main_nms. There are commentaries there and I hope that you will understand what is going on. Just a few comments.

In Fortran the first 5 symbols in each line are reserved for a numerical label, the sixth symbol, if present, shows that the previous line continues over. Then for example go to 35 means that there is jump to the label 35, do 20 i = kp1,n means that the loop goes until the line with label 20 and so on.

The array indices start from 1. This is the difference with C and C++. Another difference is that arrays in Fortran are stored column-wise.  Finally the memory allocation in Fortran is static, that is, the array dimension must be specified at the compile time. The latter leads to the next practice. In the main program one chooses the maximum dimension of the array (LDA in main_nms) and then at the run-time one can use any matrix dimension from one to LDA. This in turn means that the matrix will be written in memory  as follows. In the example LDA=10,  and  the compiler  at compile time allocates for the matrix A 100 double precision numbers continuously (just imagine one-dimensional array of 100 doubles). The first column of the matrix 3×3 employed in the example takes the first three places in this array. The second column (also three numbers) starts however from the eleventh element of the one-dimension array, and so on. It looks awkward nowadays but that’s it.

The subroutines allocate no memory by themselves. I like this practice as it makes the memory management quite clean. When a subroutine in a Fortran library needs a temporary storage, it expects to get some array for that from the main program. When a matrix is passed to a subroutine (internally this is just one-dimensional array), it is necessary to specify the dimension of the matrix N as well as LDA (the maximum number of rows). Otherwise, the subroutine will not be able to locate correctly the element A(i, j).

Now let us compile and run the code. With gfortran the command looks as follows

$ gfortran main_nms.f decomp.f solve.f -o main_nms

and this should be similar for other Fortran compilers. When you run the binary, you will see that the results should be very close with the expected solution.

In the command above gfortran makes several steps at once. It compiles files separately and then links the object files together. The disadvantage here is that when one makes changes in one file, the command recompiles all three files. It could be more efficient to split the job to pieces

$ gfortran -c decomp.f
$ gfortran -c solve.f
$ gfortran -c main_nms.f
$ gfortran main_nms.o solve.o decomp.o -o main_nms

There are more commands, but then when one changes something in the main program it will be necessary to recompile it only. One can automate the process above with a makefile. Then the command

$ make -f main_nms.make

will do the job for all the case. Try it, then make some changes in the main program and try it again.

In the makefile one also finds optimization flags. It is very important to turn optimization on, but for this simple example this actually does not matter.

Previous:

Overview
Linear Solve in Python (NumPy and SciPy)

Next:

Using decomp and solve from C++
Using LAPACK from C++