(LU, pivots) = Matrices.LU(A); (LU, pivots, info) = Matrices.LU(A);
This function call returns the LU decomposition of a "Real[m,n]" matrix A, i.e.,
P*L*U = A
where P is a permutation matrix (implicitely
defined by vector pivots
),
L is a lower triangular matrix with unit
diagonal elements (lower trapezoidal if m > n), and
U is an upper triangular matrix (upper trapezoidal if m < n).
Matrices L and U are stored in the returned
matrix LU
(the diagonal of L is not stored).
With the companion function
Matrices.LU_solve,
this decomposition can be used to solve
linear systems (P*L*U)*x = b with different right
hand side vectors b. If a linear system of equations with
just one right hand side vector b shall be solved, it is
more convenient to just use the function
Matrices.solve.
The optional third (Integer) output argument has the following meaning:
info = 0: Successful exit info > 0: If info = i then U[i,i] is exactly zero. The factorization has been completed, but the factor U is exactly singular, and division by zero will occur if it is used to solve a system of equations.
The LU factorization is computed with the LAPACK function "dgetrf", i.e., by Gaussian elemination using partial pivoting with row interchanges. Vector "pivots" are the pivot indices, i.e., for 1 ≤ i ≤ min(m,n), row i of matrix A was interchanged with row pivots[i].
Real A[3,3] = [1,2,3; 3,4,5; 2,1,4]; Real b1[3] = {10,22,12}; Real b2[3] = { 7,13,10}; Real LU[3,3]; Integer pivots[3]; Real x1[3]; Real x2[3]; algorithm (LU, pivots) := Matrices.LU(A); x1 := Matrices.LU_solve(LU, pivots, b1); // x1 = {3,2,1} x2 := Matrices.LU_solve(LU, pivots, b2); // x2 = {1,0,2}
Matrices.LU_solve, Matrices.solve
function LU extends Modelica.Icons.Function; input Real A[:, :] "Square or rectangular matrix"; output Real LU[size(A, 1), size(A, 2)] = A "L,U factors (used with LU_solve(..))"; output Integer pivots[min(size(A, 1), size(A, 2))] "Pivot indices (used with LU_solve(..))"; output Integer info "Information"; end LU;