symmetric_tridiagonal.rs - source
Description: Source of the Rust file `/Users/erlendbasso/.cargo/registry/src/github.com-1ecc6299db9ec823/nalgebra-0.32.1/src/linalg/symmetric_tridiagonal.rs`.
Keywords: rust, rustlang, rust-lang

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
    100
    101
    102
    103
    104
    105
    106
    107
    108
    109
    110
    111
    112
    113
    114
    115
    116
    117
    118
    119
    120
    121
    122
    123
    124
    125
    126
    127
    128
    129
    130
    131
    132
    133
    134
    135
    136
    137
    138
    139
    140
    141
    142
    143
    144
    145
    146
    147
    148
    149
    150
    151
    152
    153
    154
    155
    156
    157
    158
    159
    160
    161
    162
    163
    164
    165
    166
    167
    168
    169
    170

    #[cfg(feature = "serde-serialize-no-std")]
    use serde::{Deserialize, Serialize};

    use crate::allocator::Allocator;
    use crate::base::{DefaultAllocator, OMatrix, OVector};
    use crate::dimension::{Const, DimDiff, DimSub, U1};
    use simba::scalar::ComplexField;

    use crate::linalg::householder;
    use crate::Matrix;
    use std::mem::MaybeUninit;

    /// Tridiagonalization of a symmetric matrix.
    #[cfg_attr(feature = "serde-serialize-no-std", derive(Serialize, Deserialize))]
    #[cfg_attr(
        feature = "serde-serialize-no-std",
        serde(bound(serialize = "DefaultAllocator: Allocator<T, D, D> +
                               Allocator<T, DimDiff<D, U1>>,
             OMatrix<T, D, D>: Serialize,
             OVector<T, DimDiff<D, U1>>: Serialize"))
    )]
    #[cfg_attr(
        feature = "serde-serialize-no-std",
        serde(bound(deserialize = "DefaultAllocator: Allocator<T, D, D> +
                               Allocator<T, DimDiff<D, U1>>,
             OMatrix<T, D, D>: Deserialize<'de>,
             OVector<T, DimDiff<D, U1>>: Deserialize<'de>"))
    )]
    #[derive(Clone, Debug)]
    pub struct SymmetricTridiagonal<T: ComplexField, D: DimSub<U1>>
    where
        DefaultAllocator: Allocator<T, D, D> + Allocator<T, DimDiff<D, U1>>,
    {
        tri: OMatrix<T, D, D>,
        off_diagonal: OVector<T, DimDiff<D, U1>>,
    }

    impl<T: ComplexField, D: DimSub<U1>> Copy for SymmetricTridiagonal<T, D>
    where
        DefaultAllocator: Allocator<T, D, D> + Allocator<T, DimDiff<D, U1>>,
        OMatrix<T, D, D>: Copy,
        OVector<T, DimDiff<D, U1>>: Copy,
    {
    }

    impl<T: ComplexField, D: DimSub<U1>> SymmetricTridiagonal<T, D>
    where
        DefaultAllocator: Allocator<T, D, D> + Allocator<T, DimDiff<D, U1>>,
    {
        /// Computes the tridiagonalization of the symmetric matrix `m`.
        ///
        /// Only the lower-triangular part (including the diagonal) of `m` is read.
        pub fn new(mut m: OMatrix<T, D, D>) -> Self {
            let dim = m.shape_generic().0;

            assert!(
                m.is_square(),
                "Unable to compute the symmetric tridiagonal decomposition of a non-square matrix."
            );
            assert!(
                dim.value() != 0,
                "Unable to compute the symmetric tridiagonal decomposition of an empty matrix."
            );

            let mut off_diagonal = Matrix::uninit(dim.sub(Const::<1>), Const::<1>);
            let mut p = Matrix::zeros_generic(dim.sub(Const::<1>), Const::<1>);

            for i in 0..dim.value() - 1 {
                let mut m = m.rows_range_mut(i + 1..);
                let (mut axis, mut m) = m.columns_range_pair_mut(i, i + 1..);

                let (norm, not_zero) = householder::reflection_axis_mut(&mut axis);
                off_diagonal[i] = MaybeUninit::new(norm);

                if not_zero {
                    let mut p = p.rows_range_mut(i..);

                    p.hegemv(crate::convert(2.0), &m, &axis, T::zero());

                    let dot = axis.dotc(&p);
                    m.hegerc(-T::one(), &p, &axis, T::one());
                    m.hegerc(-T::one(), &axis, &p, T::one());
                    m.hegerc(dot * crate::convert(2.0), &axis, &axis, T::one());
                }
            }

            // Safety: off_diagonal has been fully initialized.
            let off_diagonal = unsafe { off_diagonal.assume_init() };
            Self {
                tri: m,
                off_diagonal,
            }
        }

        #[doc(hidden)]
        // For debugging.
        pub fn internal_tri(&self) -> &OMatrix<T, D, D> {
            &self.tri
        }

        /// Retrieve the orthogonal transformation, diagonal, and off diagonal elements of this
        /// decomposition.
        pub fn unpack(
            self,
        ) -> (
            OMatrix<T, D, D>,
            OVector<T::RealField, D>,
            OVector<T::RealField, DimDiff<D, U1>>,
        )
        where
            DefaultAllocator: Allocator<T::RealField, D> + Allocator<T::RealField, DimDiff<D, U1>>,
        {
            let diag = self.diagonal();
            let q = self.q();

            (q, diag, self.off_diagonal.map(T::modulus))
        }

        /// Retrieve the diagonal, and off diagonal elements of this decomposition.
        pub fn unpack_tridiagonal(
            self,
        ) -> (
            OVector<T::RealField, D>,
            OVector<T::RealField, DimDiff<D, U1>>,
        )
        where
            DefaultAllocator: Allocator<T::RealField, D> + Allocator<T::RealField, DimDiff<D, U1>>,
        {
            (self.diagonal(), self.off_diagonal.map(T::modulus))
        }

        /// The diagonal components of this decomposition.
        #[must_use]
        pub fn diagonal(&self) -> OVector<T::RealField, D>
        where
            DefaultAllocator: Allocator<T::RealField, D>,
        {
            self.tri.map_diagonal(|e| e.real())
        }

        /// The off-diagonal components of this decomposition.
        #[must_use]
        pub fn off_diagonal(&self) -> OVector<T::RealField, DimDiff<D, U1>>
        where
            DefaultAllocator: Allocator<T::RealField, DimDiff<D, U1>>,
        {
            self.off_diagonal.map(T::modulus)
        }

        /// Computes the orthogonal matrix `Q` of this decomposition.
        #[must_use]
        pub fn q(&self) -> OMatrix<T, D, D> {
            householder::assemble_q(&self.tri, self.off_diagonal.as_slice())
        }

        /// Recomputes the original symmetric matrix.
        pub fn recompose(mut self) -> OMatrix<T, D, D> {
            let q = self.q();
            self.tri.fill_lower_triangle(T::zero(), 2);
            self.tri.fill_upper_triangle(T::zero(), 2);

            for i in 0..self.off_diagonal.len() {
                let val = T::from_real(self.off_diagonal[i].clone().modulus());
                self.tri[(i + 1, i)] = val.clone();
                self.tri[(i, i + 1)] = val;
            }

            &q * self.tri * q.adjoint()
        }
    }