nalgebra/linalg/
balancing.rs

1//! Functions for balancing a matrix.
2
3use simba::scalar::RealField;
4use std::ops::{DivAssign, MulAssign};
5
6use crate::allocator::Allocator;
7use crate::base::dimension::Dim;
8use crate::base::{Const, DefaultAllocator, OMatrix, OVector};
9
10/// Applies in-place a modified Parlett and Reinsch matrix balancing with 2-norm to the matrix and returns
11/// the corresponding diagonal transformation.
12///
13/// See <https://arxiv.org/pdf/1401.5766.pdf>
14pub fn balance_parlett_reinsch<T: RealField, D: Dim>(matrix: &mut OMatrix<T, D, D>) -> OVector<T, D>
15where
16    DefaultAllocator: Allocator<D, D> + Allocator<D>,
17{
18    assert!(matrix.is_square(), "Unable to balance a non-square matrix.");
19
20    let dim = matrix.shape_generic().0;
21    let radix: T = crate::convert(2.0f64);
22    let mut d = OVector::from_element_generic(dim, Const::<1>, T::one());
23
24    let mut converged = false;
25
26    while !converged {
27        converged = true;
28
29        for i in 0..dim.value() {
30            let mut n_col = matrix.column(i).norm_squared();
31            let mut n_row = matrix.row(i).norm_squared();
32            let mut f = T::one();
33
34            let s = n_col.clone() + n_row.clone();
35            n_col = n_col.sqrt();
36            n_row = n_row.sqrt();
37
38            if n_col.clone().is_zero() || n_row.clone().is_zero() {
39                continue;
40            }
41
42            while n_col.clone() < n_row.clone() / radix.clone() {
43                n_col *= radix.clone();
44                n_row /= radix.clone();
45                f *= radix.clone();
46            }
47
48            while n_col.clone() >= n_row.clone() * radix.clone() {
49                n_col /= radix.clone();
50                n_row *= radix.clone();
51                f /= radix.clone();
52            }
53
54            let eps: T = crate::convert(0.95);
55            #[allow(clippy::suspicious_operation_groupings)]
56            if n_col.clone() * n_col + n_row.clone() * n_row < eps * s {
57                converged = false;
58                d[i] *= f.clone();
59                matrix.column_mut(i).mul_assign(f.clone());
60                matrix.row_mut(i).div_assign(f.clone());
61            }
62        }
63    }
64
65    d
66}
67
68/// Computes in-place `D * m * D.inverse()`, where `D` is the matrix with diagonal `d`.
69pub fn unbalance<T: RealField, D: Dim>(m: &mut OMatrix<T, D, D>, d: &OVector<T, D>)
70where
71    DefaultAllocator: Allocator<D, D> + Allocator<D>,
72{
73    assert!(m.is_square(), "Unable to unbalance a non-square matrix.");
74    assert_eq!(m.nrows(), d.len(), "Unbalancing: mismatched dimensions.");
75
76    for j in 0..d.len() {
77        let mut col = m.column_mut(j);
78        let denom = T::one() / d[j].clone();
79
80        for i in 0..d.len() {
81            col[i] *= d[i].clone() * denom.clone();
82        }
83    }
84}
nalgebra/linalg/balancing.rs

nalgebra/linalg/
balancing.rs
