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)