12.1 LU Decomposition

A general square matrix A has an LU decomposition into upper and lower triangular matrices,

P A = L U

where P is a permutation matrix, L is unit lower triangular matrix and U is upper triangular matrix. For square matrices this decomposition can be used to convert the linear system A x = b into a pair of triangular systems (L y = P b, U x = y), which can be solved by forward and back-substitution.

Function: int gsl_linalg_LU_decomp (gsl_matrix * A, gsl_permutation * p, int *signum)

Function: int gsl_linalg_complex_LU_decomp (gsl_matrix_complex * A, gsl_permutation * p, int *signum)

These functions factorize the square matrix A into the LU decomposition PA = LU. On output the diagonal and upper triangular part of the input matrix A contain the matrix U. The lower triangular part of the input matrix (excluding the diagonal) contains L. The diagonal elements of L are unity, and are not stored.

The permutation matrix P is encoded in the permutation p. The j-th column of the matrix P is given by the k-th column of the identity matrix, where k = p_j the j-th element of the permutation vector. The sign of the permutation is given by signum. It has the value (-1)^n, where n is the number of interchanges in the permutation.

The algorithm used in the decomposition is Gaussian Elimination with partial pivoting (Golub & Van Loan, Matrix Computations, Algorithm 3.4.1).

Function: int gsl_linalg_LU_solve (const gsl_matrix * LU, const gsl_permutation * p, const gsl_vector * b, gsl_vector * x)

Function: int gsl_linalg_complex_LU_solve (const gsl_matrix_complex * LU, const gsl_permutation * p, const gsl_vector_complex * b, gsl_vector_complex * x)

These functions solve the system A x = b using the LU decomposition of A into (LU, p) given by gsl_linalg_LU_decomp or gsl_linalg_complex_LU_decomp.

Function: int gsl_linalg_LU_svx (const gsl_matrix * LU, const gsl_permutation * p, gsl_vector * x)

Function: int gsl_linalg_complex_LU_svx (const gsl_matrix_complex * LU, const gsl_permutation * p, gsl_vector_complex * x)

These functions solve the system A x = b in-place using the LU decomposition of A into (LU,p). On input x should contain the right-hand side b, which is replaced by the solution on output.

Function: int gsl_linalg_LU_refine (const gsl_matrix * A, const gsl_matrix * LU, const gsl_permutation * p, const gsl_vector * b, gsl_vector * x, gsl_vector * residual)

Function: int gsl_linalg_complex_LU_refine (const gsl_matrix_complex * A, const gsl_matrix_complex * LU, const gsl_permutation * p, const gsl_vector_complex * b, gsl_vector_complex * x, gsl_vector_complex * residual)

These functions apply an iterative improvement to x, the solution of A x = b, using the LU decomposition of A into (LU,p). The initial residual r = A x - b is also computed and stored in residual.

Function: int gsl_linalg_LU_invert (const gsl_matrix * LU, const gsl_permutation * p, gsl_matrix * inverse)

Function: int gsl_complex_linalg_LU_invert (const gsl_matrix_complex * LU, const gsl_permutation * p, gsl_matrix_complex * inverse)

These functions compute the inverse of a matrix A from its LU decomposition (LU,p), storing the result in the matrix inverse. The inverse is computed by solving the system A x = b for each column of the identity matrix.

Function: double gsl_linalg_LU_det (gsl_matrix * LU, int signum)

Function: gsl_complex gsl_linalg_complex_LU_det (gsl_matrix_complex * LU, int signum)

These functions compute the determinant of a matrix A from its LU decomposition, LU. The determinant is computed as the product of the diagonal elements of U and the sign of the row permutation signum.

Function: double gsl_linalg_LU_lndet (gsl_matrix * LU)

Function: double gsl_linalg_complex_LU_lndet (gsl_matrix_complex * LU)

These functions compute the logarithm of the absolute value of the determinant of a matrix A, \ln|det(A)|, from its LU decomposition, LU. This function may be useful if the direct computation of the determinant would overflow or underflow.

Function: int gsl_linalg_LU_sgndet (gsl_matrix * LU, int signum)

Function: gsl_complex gsl_linalg_complex_LU_sgndet (gsl_matrix_complex * LU, int signum)

These functions compute the sign or phase factor of the determinant of a matrix A, det(A)/|det(A)|, from its LU decomposition, LU.