LU decomposition: Difference between revisions
Content added Content deleted
(Added ATS.) |
No edit summary |
||
Line 5,312: | Line 5,312: | ||
[0, 0, 0, 1]] |
[0, 0, 0, 1]] |
||
</pre> |
</pre> |
||
===Alternative with abstalg==== |
|||
{{libheader| abstalg}} |
|||
This one implements a naive LU decomposition with arbitrary fields, so it works over Rational Numbers as well as floats, (or any field) |
|||
<syntaxhighlight lang="rust"> |
|||
use abstalg::*; |
|||
pub struct Matrix2D<'a, F> |
|||
where |
|||
F: Field, |
|||
{ |
|||
field: MatrixRing<F>, |
|||
data: &'a mut Vec<<F as Domain>::Elem>, |
|||
rows: usize, |
|||
cols: usize, |
|||
} |
|||
impl<'a, F: Field + Clone> Matrix2D<'a, F> { |
|||
pub fn new(field: F, data: &'a mut Vec<<F as Domain>::Elem>, rows: usize, cols: usize) -> Self { |
|||
assert_eq!(rows * cols, data.len(), "Data does not match dimensions"); |
|||
Matrix2D { |
|||
field: MatrixRing::<F>::new(field, rows), |
|||
data, |
|||
rows, |
|||
cols, |
|||
} |
|||
} |
|||
pub fn get(&self, row: usize, col: usize) -> &<F as Domain>::Elem { |
|||
assert!(row < self.rows && col < self.cols, "Index out of bounds"); |
|||
&self.data[row * self.cols + col] |
|||
} |
|||
pub fn get_mut(&mut self, row: usize, col: usize) -> &mut <F as Domain>::Elem { |
|||
assert!(row < self.rows && col < self.cols, "Index out of bounds"); |
|||
&mut self.data[row * self.cols + col] |
|||
} |
|||
pub fn get_row(&self, row: usize) -> Vec<<F as Domain>::Elem> { |
|||
assert!(row < self.rows, "Row index out of bounds"); |
|||
let mut result = Vec::new(); |
|||
for col in 0..self.cols { |
|||
result.push(self.get(row, col).clone()); |
|||
} |
|||
result |
|||
} |
|||
pub fn get_col(&self, col: usize) -> Vec<<F as Domain>::Elem> { |
|||
assert!(col < self.cols, "Column index out of bounds"); |
|||
let mut result = Vec::new(); |
|||
for row in 0..self.rows { |
|||
result.push(self.get(row, col).clone()); |
|||
} |
|||
result |
|||
} |
|||
pub fn set_row(&mut self, row: usize, new_row: Vec<<F as Domain>::Elem>) { |
|||
assert!(row < self.rows, "Row index out of bounds"); |
|||
assert_eq!(new_row.len(), self.cols, "New row has wrong length"); |
|||
for col in 0..self.cols { |
|||
*self.get_mut(row, col) = new_row[col].clone(); |
|||
} |
|||
} |
|||
pub fn set_col(&mut self, col: usize, new_col: Vec<<F as Domain>::Elem>) { |
|||
assert!(col < self.cols, "Column index out of bounds"); |
|||
assert_eq!(new_col.len(), self.rows, "New column has wrong length"); |
|||
for row in 0..self.rows { |
|||
*self.get_mut(row, col) = new_col[row].clone(); |
|||
} |
|||
} |
|||
pub fn swap_rows(&mut self, row1: usize, row2: usize) { |
|||
assert!( |
|||
row1 < self.rows && row2 < self.rows, |
|||
"Row index out of bounds" |
|||
); |
|||
if row1 != row2 { |
|||
for col in 0..self.cols { |
|||
let temp = self.get(row1, col).clone(); |
|||
*self.get_mut(row1, col) = self.get(row2, col).clone(); |
|||
*self.get_mut(row2, col) = temp; |
|||
} |
|||
} |
|||
} |
|||
pub fn l_u_decomposition(&mut self) -> Result<Vec<<F as Domain>::Elem>, String> |
|||
where |
|||
F: Clone, |
|||
{ |
|||
// Let base = field.base() |
|||
let base = self.field.base().clone(); |
|||
// Let v_a = VectorAlgebra(base, cols) |
|||
let v_a = VectorAlgebra::new(base.clone(), self.cols); |
|||
// Let the_l_matrix = I (creates an identity matrix) |
|||
let mut the_l_matrix: Vec<_> = self.field.int(1); |
|||
// Let l_matrix = Matrix2D(base, the_l_matrix, rows, cols) |
|||
let mut l_matrix = Matrix2D::new(base.clone(), &mut the_l_matrix, self.rows, self.cols); |
|||
// For each pivot in min(rows, cols) |
|||
for pivot in 0..std::cmp::min(self.rows, self.cols) { |
|||
// Let pivot_row = self.get_row(pivot) |
|||
let pivot_row = self.get_row(pivot); |
|||
// If pivot element (pivot_row[pivot]) is zero, LU decomposition is not possible |
|||
if base.is_zero(&pivot_row[pivot]) { |
|||
return Err( |
|||
"LU decomposition without pivoting is not possible for this matrix".into(), |
|||
); |
|||
} |
|||
// Let pivot_entry_inv = 1 / pivot_row[pivot] |
|||
let pivot_entry_inv = base.inv(&pivot_row[pivot]); |
|||
// For each row_idx in (pivot + 1) to rows |
|||
for row_idx in (pivot + 1)..self.rows { |
|||
// Let row = self.get_row(row_idx) |
|||
let mut row = self.get_row(row_idx); |
|||
// Let scale = row[pivot] * pivot_entry_inv |
|||
let scale = base.mul(&row[pivot], &pivot_entry_inv); |
|||
// row += -scale * pivot_row (Vector addition and scalar multiplication) |
|||
v_a.add_assign( |
|||
&mut row, |
|||
&v_a.neg(&mul_vector(&v_a, scale.clone(), pivot_row.clone())), |
|||
); |
|||
// l_matrix[row_idx][pivot] = scale |
|||
*l_matrix.get_mut(row_idx, pivot) = scale; |
|||
// self.set_row(row_idx, row) (Sets the modified row back into the matrix) |
|||
self.set_row(row_idx, row); |
|||
} |
|||
} |
|||
// Returns the L matrix |
|||
Ok(the_l_matrix) |
|||
} |
|||
pub fn p_l_u_decomposition( |
|||
&self, |
|||
) -> Result< |
|||
( |
|||
Vec<<F as Domain>::Elem>, |
|||
Vec<<F as Domain>::Elem>, |
|||
Vec<<F as Domain>::Elem>, |
|||
), |
|||
String, |
|||
> |
|||
where |
|||
F: Clone, |
|||
{ |
|||
let base = self.field.base().clone(); |
|||
let mut self2 = (*self.data).clone(); |
|||
let mut cloned_vector = Matrix2D::new(base.clone(), &mut self2, self.rows, self.cols); |
|||
let mut pivot_row = 0; |
|||
let mut the_p_matrix: Vec<_> = self.field.zero(); |
|||
let mut p_matrix = Matrix2D::new(base.clone(), &mut the_p_matrix, self.rows, self.cols); |
|||
//let mut u_matrix = self.clone(); //Initializes the U matrix as a copy of the original matrix |
|||
for pivot_col in 0..self.cols { |
|||
// Find a non-zero entry in the pivot column |
|||
let swap_row = (pivot_row..self.rows) |
|||
.find(|&row| !base.equals(cloned_vector.get(row, pivot_col), &base.zero())); |
|||
match swap_row { |
|||
Some(swap_row) => { |
|||
// Swap rows in U and P matrices to bring the non-zero entry to the pivot position |
|||
cloned_vector.swap_rows(pivot_row, swap_row); |
|||
p_matrix.swap_rows(pivot_row, swap_row); |
|||
pivot_row += 1; |
|||
} |
|||
None => { |
|||
// If there are no non-zero entries in the pivot column, just proceed to the next column |
|||
} |
|||
} |
|||
} |
|||
// Set the diagonals of P to 1 |
|||
for i in 0..self.rows { |
|||
*p_matrix.get_mut(i, i) = base.one(); |
|||
} |
|||
// Run the LU decomposition on the permuted U matrix |
|||
let l_u_result = cloned_vector.l_u_decomposition(); |
|||
match l_u_result { |
|||
Ok(the_l_matrix) => Ok((the_p_matrix, the_l_matrix, cloned_vector.data.clone())), |
|||
Err(e) => Err(e), |
|||
} |
|||
} |
|||
} |
|||
use std::{error::Error, fmt}; |
|||
impl<'a, T> fmt::Display for Matrix2D<'a, T> |
|||
where |
|||
<T as Domain>::Elem: fmt::Display, |
|||
T: Field, |
|||
{ |
|||
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { |
|||
for row in 0..self.rows { |
|||
for col in 0..self.cols { |
|||
write!(f, "{} ", self.get(row, col))?; |
|||
} |
|||
writeln!(f)?; |
|||
} |
|||
Ok(()) |
|||
} |
|||
} |
|||
fn mul_vector<T>( |
|||
f: &VectorAlgebra<T>, |
|||
a: <T as Domain>::Elem, |
|||
d: Vec<<T as Domain>::Elem>, |
|||
) -> <VectorAlgebra<T> as Domain>::Elem |
|||
where |
|||
T: Field, |
|||
{ |
|||
f.mul(&d, &f.diagonal(a)) |
|||
} |
|||
#[cfg(test)] |
|||
mod tests { |
|||
use super::*; // bring into scope everything from the parent module |
|||
use abstalg::ReducedFractions; |
|||
use abstalg::I32; |
|||
fn test_p_l_u_decomposition(matrix: Vec<isize>, size: usize) { |
|||
// This test assumes that the base field is rational numbers. |
|||
// Create a test 4x4 matrix |
|||
let (matrix, size): (Vec<isize>, usize) = |
|||
(vec![11, 9, 24, 2, 1, 5, 2, 6, 3, 17, 18, 1, 2, 5, 7, 1], 4); |
|||
let field = abstalg::ReducedFractions::new(abstalg::I32); |
|||
let matrix_ring = MatrixRing::new(field.clone(), size); |
|||
let mut matrix: Vec<_> = matrix.clone().into_iter().map(|i| field.int(i)).collect(); |
|||
let mut matrix2d = Matrix2D::new(field.clone(), &mut matrix, size, size); |
|||
// Decompose the matrix using the p_l_u_decomposition function |
|||
let p_l_u_decomposition_result = matrix2d.p_l_u_decomposition().unwrap(); |
|||
let (mut p_matrix, mut l_matrix, mut u_matrix) = p_l_u_decomposition_result; |
|||
// Convert the matrices back to Matrix2D form for printing |
|||
let p_matrix2d = Matrix2D::new(field.clone(), &mut p_matrix, size, size); |
|||
let l_matrix2d = Matrix2D::new(field.clone(), &mut l_matrix, size, size); |
|||
let u_matrix2d = Matrix2D::new(field.clone(), &mut u_matrix, size, size); |
|||
println!("P={} L={} U={}", p_matrix2d, l_matrix2d, u_matrix2d,); |
|||
// Multiply the resulting P, L, and U matrices |
|||
let p_l = matrix_ring.mul(&p_matrix, &l_matrix); |
|||
let mut p_l_u = matrix_ring.mul(&p_l, &u_matrix); |
|||
//let p_l_u_2d = Matrix2D::new(field.clone(), &mut p_l_u, 4, 4); |
|||
// Check that the product of P, L, and U is equal to the original matrix |
|||
assert_eq!(matrix, p_l_u); |
|||
} |
|||
#[test] |
|||
fn test_p_l_u_decomposition_example() { |
|||
test_p_l_u_decomposition(vec![ |
|||
11, 9, 24, 2, |
|||
1, 5, 2, 6, |
|||
3, 17, 18, 1, |
|||
2, 5, 7, 1, |
|||
], 4); |
|||
test_p_l_u_decomposition(vec![ |
|||
1, 3, 5, |
|||
2, 4, 7, |
|||
1, 1, 0, |
|||
], 3); |
|||
} |
|||
} |
|||
</syntaxhighlight> |
|||
=={{header|Sidef}}== |
=={{header|Sidef}}== |
||
{{trans|Raku}} |
{{trans|Raku}} |