pub fn LU_decomp<T>(A: &Array2<T>) -> (Array2<T>, Array2<T>, Array2<T>)
where
    T: NdFloat,
{
    let matrix_dimension = A.nrows();
    assert_eq!(
        matrix_dimension,
        A.ncols(),
        "Tried LU decomposition with a non-square matrix."
    );
    let P = pivot(&A);
    let pivotized_A = P.dot(A);

    let mut L: Array2<T> = Array::eye(matrix_dimension);
    let mut U: Array2<T> = Array::zeros((matrix_dimension, matrix_dimension));
    for idx_col in 0..matrix_dimension {
        // fill U
        for idx_row in 0..idx_col + 1 {
            U[[idx_row, idx_col]] = pivotized_A[[idx_row, idx_col]]
                - U.slice(s![0..idx_row, idx_col])
                    .dot(&L.slice(s![idx_row, 0..idx_row]));
        }
        // fill L
        for idx_row in idx_col + 1..matrix_dimension {
            L[[idx_row, idx_col]] = (pivotized_A[[idx_row, idx_col]]
                - U.slice(s![0..idx_col, idx_col])
                    .dot(&L.slice(s![idx_row, 0..idx_col])))
                / U[[idx_col, idx_col]];
        }
    }
    (L, U, P)
}
/// Performs partial pivoted LU decomposition of A such that P A = LU
/// with L a lower triangular matrix and U an upper triangular matrix
/// A needs to be a square matrix
Performs partial pivoted LU decomposition of A such that P A = LU
with L a lower triangular matrix and U an upper triangular matrix
A needs to be a square matrix
Performs partial pivoted LU decomposition of A such that P A = LU with L a lower triangular matrix and U an upper triangular matrix A needs to be a square matrix
pub fn LU_decomp<T>(A: &Array2<T>) -> (Array2<T>, Array2<T>, Array2<T>)
where
    T: NdFloat,
{
    let matrix_dimension = A.nrows();
    assert_eq!(
        matrix_dimension,
        A.ncols(),
        "Tried LU decomposition with a non-square matrix."
    );
    let P = pivot(&A);
    let pivotized_A = P.dot(A);
    let mut L: Array2<T> = Array::eye(matrix_dimension);
    let mut U: Array2<T> = Array::zeros((matrix_dimension, matrix_dimension));
    for idx_col in 0..matrix_dimension {
        // fill U
        for idx_row in 0..idx_col + 1 {
            U[[idx_row, idx_col]] = pivotized_A[[idx_row, idx_col]]
                - U.slice(s![0..idx_row, idx_col])
                    .dot(&L.slice(s![idx_row, 0..idx_row]));
        }
        // fill L
        for idx_row in idx_col + 1..matrix_dimension {
            L[[idx_row, idx_col]] = (pivotized_A[[idx_row, idx_col]]
                - U.slice(s![0..idx_col, idx_col])
                    .dot(&L.slice(s![idx_row, 0..idx_col])))
                / U[[idx_col, idx_col]];
        }
    }
    (L, U, P)
}
pub fn LU_decomp < T > ( A : & Array2 < T > ) -> ( Array2 < T > , Array2 < T > , Array2 < T > ) where T : NdFloat , { let matrix_dimension = A . nrows ( ) ; assert_eq ! ( matrix_dimension A ncols ( ) " " ) ; let P = pivot ( & A ) ; let pivotized_A = P . dot ( A ) ; let mut L : Array2 < T > = Array :: eye ( matrix_dimension ) ; let mut U : Array2 < T > = Array :: zeros ( ( matrix_dimension , matrix_dimension ) ) ; for idx_col in 0 .. matrix_dimension { // fill U for idx_row in 0 .. idx_col + 1 { U [ [ idx_row , idx_col ] ] = pivotized_A [ [ idx_row , idx_col ] ] - U . slice ( s ! [ 0 idx_row idx_col ] ) . dot ( & L . slice ( s ! [ idx_row 0 idx_row ] ) ) ; } // fill L for idx_row in idx_col + 1 .. matrix_dimension { L [ [ idx_row , idx_col ] ] = ( pivotized_A [ [ idx_row , idx_col ] ] - U . slice ( s ! [ 0 idx_col idx_col ] ) . dot ( & L . slice ( s ! [ idx_row 0 idx_col ] ) ) ) / U [ [ idx_col , idx_col ] ] ; } } ( L , U , P ) }
Performs partial pivoted LU decomposition of A such that P A = LU
with L a lower triangular matrix and U an upper triangular matrix
A needs to be a square matrix
Performs partial pivoted LU decomposition of A such that P A = LU with L a lower triangular matrix and U an upper triangular matrix A needs to be a square matrix
A
&Array2<T>
A
&Array2<T>