LU Decomposition

A general square matrix @math{A} has an @math{LU} decomposition into upper and lower triangular matrices,

where @math{P} is a permutation matrix, @math{L} is unit lower triangular matrix and @math{U} is upper triangular matrix. For square matrices this decomposition can be used to convert the linear system @math{A x = b} into a pair of triangular systems (@math{L y = P b}, @math{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 @math{LU} decomposition @math{PA = LU}. On output the diagonal and upper triangular part of the input matrix A contain the matrix @math{U}. The lower triangular part of the input matrix (excluding the diagonal) contains @math{L}. The diagonal elements of @math{L} are unity, and are not stored.

The permutation matrix @math{P} is encoded in the permutation p. The @math{j}-th column of the matrix @math{P} is given by the @math{k}-th column of the identity matrix, where @math{k = p_j} the @math{j}-th element of the permutation vector. The sign of the permutation is given by signum. It has the value @math{(-1)^n}, where @math{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 @math{A x = b} using the @math{LU} decomposition of @math{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 @math{A x = b} in-place using the @math{LU} decomposition of @math{A} into (LU,p). On input x should contain the right-hand side @math{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 @math{A x = b}, using the @math{LU} decomposition of @math{A} into (LU,p). The initial residual @math{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 @math{A} from its @math{LU} decomposition (LU,p), storing the result in the matrix inverse. The inverse is computed by solving the system @math{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 @math{A} from its @math{LU} decomposition, LU. The determinant is computed as the product of the diagonal elements of @math{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 @math{A}, @math{\ln|det(A)|}, from its @math{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 @math{A}, @math{det(A)/|det(A)|}, from its @math{LU} decomposition, LU.