Solution of a system of linear equations is the task encountered most often in linear algebra. A good starting point from the theoretical side is the page in Wikipedia

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

Numerical methods to solve a linear system are divided in two groups: direct and iterative methods. For dense matrices, people use mostly direct methods and I will limit myself to them only.  The direct methods are based on the Gaussian elimination. For an unsymmetric matrix, this is the LU decomposition, for a symmetric positive definite matrix is the Cholesky decomposition, for an indefinite symmetric matrix is the LDL^T decomposition. No doubt, one could use the LU decomposition for all matrices but then it takes some more time.

The solution of a linear system consists of the two steps, first the decomposition, second the back substitution. The first step takes the most time and when it is necessary to solve a linear system with many right hand sides, one could save a lot of time by making the decomposition once and then running many back substitutions.

The most famous numerical library is LAPACK and our goal is to see how to use it in practice. We start with Python where the LAPACK library is interfaced through SciPy. This is the most comfortable way to deal with the linear algebra as one needs nothing to compile. Yet, if you develop your own system, it might be good to interface LAPACK directly.

LAPACK is a Fortran library. To learn how to interface it, I begin with the simpler code (decomp.f and solve.f) from the Forsythe book. The advantage here is that the code is quite simple and everything will be much easier. Then I show how to interface decomp.f and solve.f from C++ and only after that I will consider interfacing LAPACK.

Next:

Linear Solve in Python (NumPy and SciPy)
Using decomp and solve from Fortran
Using decomp and solve from C++
Using LAPACK from C++