AXB - Functions to solve a linear system of equations Ax=b by LU
	decomposition, invert a square matrix or directly multiply an
	inverse matrix by another matrix (without explicitely computing
	the inverse).

LU_decomposition	Decompose a matrix (A) into a lower triangular (L)
			and an upper triangular (U) such that A=LU

backward_substitution	Apply backward substitution to an LU decomposed
			matrix to solve the linear system of equations Ax=b

inverse_matrix		compute the inverse of a square non-singular matrix

inverse_matrix_multiply	computes the product A^(-1)*B without explicitely
			computing the inverse matrix

Function prototypes:
void LU_decomposition (int nrows, float **matrix, int *idx, float *d);
void backward_substitution (int nrows, float **matrix, int *idx, float *b);
void inverse_matrix (int nrows, float **matrix);
void inverse_matrix_multiply (int nrows1, float **matrix1, int ncols2,
        int nrows2, float **matrix2, float **out_matrix);

LU_decomposition:
Input:
nrows		number of rows of matrix to invert
matrix		matrix of coefficients in linear system Ax=b 

Output:
matrix		matrix containing LU decomposition (original matrix destroyed)
idx		vector recording the row permutations effected by partial
		pivoting
d		+/- 1 depending on whether the number of row interchanges
		was even or odd
backward_substitution
Input:
nrows		number of rows (and columns) of input matrix
matrix		matrix of coefficients (after LU decomposition)
idx		permutation vector obtained from routine LU_decomposition 
b		right hand side vector in equation Ax=b

Output:
b		vector with the solution
inverse_matrix
Input:
nrows		number of rows (and columns) of input matrix
matrix		matrix to invert

Output:
matrix		inverse of input matrix 
inverse_matrix_multiply
nrows1          number of rows (and columns) of matrix to invert
matrix1         square matrix to invert
ncols2          number of coulmns of second matrix
nrows2		number of rows of second matrix
matrix          second matrix (multiplicator)

Output Parameters:
out_matrix      matrix containing the product of the inverse of the first
                matrix by the second one.
Note:
matrix1 and matrix2 are not destroyed, (not clobbered)
Notes:
To solve the set of linear equations Ax=b, first do the LU decomposition of
A (which will clobber A with its LU decomposition) and then do the backward 
substitution with this new matrix and the right-hand side vector b. The vector
b will be clobbered with the solution. Both, the original matrix and vector B,
will have been destroyed.

The LU decomposition is carried out with the Crout's method with implicit
partial pivoting that guaratees that the maximum pivot is used in every
step of the algorithm.

The operation count to solve a linear system of equations via LU decomposition
is 1/3N^3 and is a factor of 3 better than the standard Gauss-Jordan algorithm
To invert a matrix the count is the same with both algorithms: N^3.

Once a linear system Ax=b has been solved, to solve another linear system
with the same matrix A but with different vetor b, ONLY the back substitution 
has to be repeated with the new b (remember that the matrix in backsubstitution
is not the original matrix but its LU decomposition)

If you want to compute A^(-1)*B from matrices A and B, it is better to 
use the subroutine inverse_matrix_multiply rather than explicitely computing
the inverse. This saves a whole martix multiplication and is also more accurate.

Refferences:
Press, Teukolsky, Vettering and Flannery, Numerical Recipes in C: 
	The art of scientific computing. Cambridge University Press.
	second edition. (1992).
Golub and Van Loan, Matrix Computations. John Hopkins University Press.
	Second Edition. (1989). 
Horn and Johnson, Matrix Analysis. Cambridge University Press. (1985).
Credits:
Adapted from discussions in Numerical Recipes, by Gabriel Alvarez (1995)