M-Labs
/
nalgebra
3
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

UDU impl: using 0-index nomenclature 

Signed-off-by: Christopher Rabotin <christopher.rabotin@gmail.com>
This commit is contained in:
Christopher Rabotin 2020-09-25 21:29:46 -06:00 committed by Crozet Sébastien
parent a8d40423ea
commit 5a7ed61e9b
1 changed files with 8 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
src/linalg/udu.rs
View File

@ -37,24 +37,24 @@ where
let mut d = MatrixN::<N, D>::zeros(); let mut d = MatrixN::<N, D>::zeros();
let mut u = MatrixN::<N, D>::zeros(); let mut u = MatrixN::<N, D>::zeros();
let n = p.ncols() - 1; let n = p.ncols();
d[(n, n)] = p[(n, n)]; d[(n - 1, n - 1)] = p[(n - 1, n - 1)];
u[(n, n)] = N::one(); u[(n - 1, n - 1)] = N::one();
for j in (0..n).rev() { for j in (0..n - 1).rev() {
u[(j, n)] = p[(j, n)] / d[(n, n)]; u[(j, n - 1)] = p[(j, n - 1)] / d[(n - 1, n - 1)];
} }
for j in (0..n).rev() { for j in (0..n - 1).rev() {
for k in j + 1..=n { for k in j + 1..n {
d[(j, j)] = d[(j, j)] + d[(k, k)] * u[(j, k)].powi(2); d[(j, j)] = d[(j, j)] + d[(k, k)] * u[(j, k)].powi(2);
} }
d[(j, j)] = p[(j, j)] - d[(j, j)]; d[(j, j)] = p[(j, j)] - d[(j, j)];
for i in (0..=j).rev() { for i in (0..=j).rev() {
for k in j + 1..=n { for k in j + 1..n {
u[(i, j)] = u[(i, j)] + d[(k, k)] * u[(j, k)] * u[(i, k)]; u[(i, j)] = u[(i, j)] + d[(k, k)] * u[(j, k)] * u[(i, k)];
} }