faer_lu/partial_pivoting/
reconstruct.rs

1use dyn_stack::{PodStack, SizeOverflow, StackReq};
2use faer_core::{
3    assert,
4    mul::triangular,
5    permutation::{Index, PermutationRef},
6    temp_mat_req, temp_mat_uninit, ComplexField, Entity, MatMut, MatRef, Parallelism,
7};
8use reborrow::*;
9use triangular::BlockStructure;
10
11#[track_caller]
12fn reconstruct_impl<I: Index, E: ComplexField>(
13    dst: MatMut<'_, E>,
14    lu_factors: Option<MatRef<'_, E>>,
15    row_perm: PermutationRef<'_, I, E>,
16    parallelism: Parallelism,
17    stack: PodStack<'_>,
18) {
19    let lu_factors = match lu_factors {
20        Some(lu_factors) => lu_factors,
21        None => dst.rb(),
22    };
23
24    let m = lu_factors.nrows();
25    let n = lu_factors.ncols();
26    let size = Ord::min(m, n);
27
28    let (mut lu, _) = temp_mat_uninit::<E>(m, n, stack);
29    let mut lu = lu.as_mut();
30
31    let (l_top, _, l_bot, _) = lu_factors.split_at(size, size);
32    let (u_left, u_right, _, _) = lu_factors.split_at(size, size);
33
34    let (lu_topleft, lu_topright, lu_botleft, _) = lu.rb_mut().split_at_mut(size, size);
35
36    triangular::matmul(
37        lu_topleft,
38        BlockStructure::Rectangular,
39        l_top,
40        BlockStructure::UnitTriangularLower,
41        u_left,
42        BlockStructure::TriangularUpper,
43        None,
44        E::faer_one(),
45        parallelism,
46    );
47    triangular::matmul(
48        lu_topright,
49        BlockStructure::Rectangular,
50        l_top,
51        BlockStructure::UnitTriangularLower,
52        u_right,
53        BlockStructure::Rectangular,
54        None,
55        E::faer_one(),
56        parallelism,
57    );
58    triangular::matmul(
59        lu_botleft,
60        BlockStructure::Rectangular,
61        l_bot,
62        BlockStructure::Rectangular,
63        u_left,
64        BlockStructure::TriangularUpper,
65        None,
66        E::faer_one(),
67        parallelism,
68    );
69
70    faer_core::permutation::permute_rows(dst, lu.rb(), row_perm.inverse());
71}
72
73/// Computes the reconstructed matrix, given its partial pivoting LU decomposition,
74/// and stores the result in `dst`.
75///
76/// # Panics
77///
78/// - Panics if the row permutation doesn't have the same dimension as the number of rows of the
79/// matrix.
80/// - Panics if the destination shape doesn't match the shape of the matrix.
81/// - Panics if the provided memory in `stack` is insufficient (see [`reconstruct_req`]).
82#[track_caller]
83pub fn reconstruct<I: Index, E: ComplexField>(
84    dst: MatMut<'_, E>,
85    lu_factors: MatRef<'_, E>,
86    row_perm: PermutationRef<'_, I, E>,
87    parallelism: Parallelism,
88    stack: PodStack<'_>,
89) {
90    assert!((dst.nrows(), dst.ncols()) == (lu_factors.nrows(), lu_factors.ncols()));
91    assert!(row_perm.len() == lu_factors.nrows());
92    reconstruct_impl(dst, Some(lu_factors), row_perm, parallelism, stack)
93}
94
95/// Computes the reconstructed matrix, given its partial pivoting LU decomposition, and stores the
96/// result in `lu_factors`.
97///
98/// # Panics
99///
100/// - Panics if the row permutation doesn't have the same dimension as the number of rows of the
101/// matrix.
102/// - Panics if the provided memory in `stack` is insufficient (see [`reconstruct_in_place_req`]).
103#[track_caller]
104pub fn reconstruct_in_place<I: Index, E: ComplexField>(
105    lu_factors: MatMut<'_, E>,
106    row_perm: PermutationRef<'_, I, E>,
107    parallelism: Parallelism,
108    stack: PodStack<'_>,
109) {
110    assert!(row_perm.len() == lu_factors.nrows());
111    reconstruct_impl(lu_factors, None, row_perm, parallelism, stack)
112}
113
114/// Computes the size and alignment of required workspace for reconstructing a matrix out of place,
115/// given its partial pivoting LU decomposition.
116pub fn reconstruct_req<I: Index, E: Entity>(
117    nrows: usize,
118    ncols: usize,
119    parallelism: Parallelism,
120) -> Result<StackReq, SizeOverflow> {
121    let _ = parallelism;
122    temp_mat_req::<E>(nrows, ncols)
123}
124
125/// Computes the size and alignment of required workspace for reconstructing a matrix in place,
126/// given its partial pivoting LU decomposition.
127pub fn reconstruct_in_place_req<I: Index, E: Entity>(
128    nrows: usize,
129    ncols: usize,
130    parallelism: Parallelism,
131) -> Result<StackReq, SizeOverflow> {
132    reconstruct_req::<I, E>(nrows, ncols, parallelism)
133}
faer_lu/partial_pivoting/reconstruct.rs

faer_lu/partial_pivoting/
reconstruct.rs
