opensrdk_linear_algebra/matrix/st/
mod.rs

1use crate::number::*;
2use crate::{matrix::ge::Matrix, matrix::*};
3use rayon::prelude::*;
4use serde::{Deserialize, Serialize};
5
6pub mod pt;
7
8pub mod ev;
9pub mod evd;
10
11#[derive(Clone, Debug, Default, PartialEq, Hash, Serialize, Deserialize)]
12pub struct SymmetricTridiagonalMatrix<T = f64>
13where
14    T: Number,
15{
16    d: Vec<T>,
17    e: Vec<T>,
18}
19
20impl<T> SymmetricTridiagonalMatrix<T>
21where
22    T: Number,
23{
24    pub fn new(dim: usize) -> Self {
25        Self {
26            d: vec![T::default(); dim],
27            e: vec![T::default(); dim.max(1) - 1],
28        }
29    }
30
31    /// - `d`: Diagonal elements. The length must be `dimension`.
32    /// - `e`: First both superdiagonal and subdiagonal elements. The length must be `dimension - 1`.
33    pub fn from(d: Vec<T>, e: Vec<T>) -> Result<Self, MatrixError> {
34        if d.len().max(1) - 1 != e.len() {
35            return Err(MatrixError::DimensionMismatch);
36        }
37
38        Ok(Self { d, e })
39    }
40
41    /// Dimension.
42    pub fn n(&self) -> usize {
43        self.d.len()
44    }
45
46    /// Diagonal elements.
47    pub fn d(&self) -> &[T] {
48        &self.d
49    }
50
51    /// First both superdiagonal and subdiagonal elements/
52    pub fn e(&self) -> &[T] {
53        &self.e
54    }
55
56    /// Returns `(self.d, self.e)`
57    pub fn eject(self) -> (Vec<T>, Vec<T>) {
58        (self.d, self.e)
59    }
60}
61
62impl SymmetricTridiagonalMatrix {
63    pub fn mat(&self) -> Matrix {
64        let n = self.d.len();
65        let mut mat = Matrix::new(n, n);
66
67        // for i in 0..n {
68        //     mat[i][i] = self.d[i];
69        // }
70
71        // for i in 0..n - 1 {
72        //     mat[i][i + 1] = self.e[i];
73        //     mat[i + 1][i] = self.e[i];
74        // }
75
76        mat.elems_mut()
77            .par_iter_mut()
78            .enumerate()
79            .map(|(k, elem)| ((k / n, k % n), elem))
80            .for_each(|((i, j), elem)| {
81                if i == j {
82                    *elem = self.d[i];
83                } else if i + 1 == j {
84                    *elem = self.e[i];
85                } else if i == j + 1 {
86                    *elem = self.e[j];
87                }
88            });
89
90        mat
91    }
92}
93
94impl SymmetricTridiagonalMatrix<c64> {
95    pub fn mat(&self, hermite: bool) -> Matrix<c64> {
96        let n = self.d.len();
97        let mut mat = Matrix::new(n, n);
98
99        mat.elems_mut()
100            .par_iter_mut()
101            .enumerate()
102            .map(|(k, elem)| ((k / n, k % n), elem))
103            .for_each(|((i, j), elem)| {
104                if i == j {
105                    *elem = self.d[i];
106                } else if i + 1 == j {
107                    *elem = self.e[i];
108                } else if i == j + 1 {
109                    *elem = if hermite { self.e[j].conj() } else { self.e[j] };
110                }
111            });
112
113        mat
114    }
115}
opensrdk_linear_algebra/matrix/st/mod.rs

opensrdk_linear_algebra/matrix/st/
mod.rs
