opensrdk_linear_algebra/matrix/bd/
mod.rs

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

opensrdk_linear_algebra/matrix/bd/
mod.rs
