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
use crate::number::*;
use crate::{matrix::ge::Matrix, matrix::*};
use rayon::prelude::*;
pub mod pt;
pub mod ev;
pub mod evd;
#[derive(Clone, Debug, Default, PartialEq, Hash)]
pub struct SymmetricTridiagonalMatrix<T = f64>
where
T: Number,
{
d: Vec<T>,
e: Vec<T>,
}
impl<T> SymmetricTridiagonalMatrix<T>
where
T: Number,
{
pub fn new(dim: usize) -> Self {
Self {
d: vec![T::default(); dim],
e: vec![T::default(); dim.max(1) - 1],
}
}
pub fn from(d: Vec<T>, e: Vec<T>) -> Result<Self, MatrixError> {
if d.len().max(1) - 1 != e.len() {
return Err(MatrixError::DimensionMismatch);
}
Ok(Self { d, e })
}
pub fn n(&self) -> usize {
self.d.len()
}
pub fn d(&self) -> &[T] {
&self.d
}
pub fn e(&self) -> &[T] {
&self.e
}
pub fn eject(self) -> (Vec<T>, Vec<T>) {
(self.d, self.e)
}
}
impl SymmetricTridiagonalMatrix {
pub fn mat(&self) -> Matrix {
let n = self.d.len();
let mut mat = Matrix::new(n, n);
mat.elems_mut()
.par_iter_mut()
.enumerate()
.map(|(k, elem)| ((k / n, k % n), elem))
.for_each(|((i, j), elem)| {
if i == j {
*elem = self.d[i];
} else if i + 1 == j {
*elem = self.e[i];
} else if i == j + 1 {
*elem = self.e[j];
}
});
mat
}
}
impl SymmetricTridiagonalMatrix<c64> {
pub fn mat(&self, hermite: bool) -> Matrix<c64> {
let n = self.d.len();
let mut mat = Matrix::new(n, n);
mat.elems_mut()
.par_iter_mut()
.enumerate()
.map(|(k, elem)| ((k / n, k % n), elem))
.for_each(|((i, j), elem)| {
if i == j {
*elem = self.d[i];
} else if i + 1 == j {
*elem = self.e[i];
} else if i == j + 1 {
*elem = if hermite { self.e[j].conj() } else { self.e[j] };
}
});
mat
}
}