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)