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