single_svdlib/
utils.rs

1use rayon::iter::ParallelIterator;
2use nalgebra_sparse::na::{DMatrix, DVector};
3use ndarray::{Array1, Array2, ShapeBuilder};
4use num_traits::{Float, Zero};
5use rayon::prelude::{IntoParallelIterator, IndexedParallelIterator};
6use single_utilities::traits::FloatOpsTS;
7use std::fmt::Debug;
8use nalgebra::{Dim, Dyn, Scalar};
9
10pub fn determine_chunk_size(nrows: usize) -> usize {
11    let num_threads = rayon::current_num_threads();
12
13    let min_rows_per_thread = 64;
14    let desired_chunks_per_thread = 4;
15
16    let target_total_chunks = num_threads * desired_chunks_per_thread;
17    let chunk_size = nrows.div_ceil(target_total_chunks);
18
19    chunk_size.max(min_rows_per_thread)
20}
21
22pub trait SMat<T: Float> {
23    fn nrows(&self) -> usize;
24    fn ncols(&self) -> usize;
25    fn nnz(&self) -> usize;
26    fn svd_opa(&self, x: &[T], y: &mut [T], transposed: bool); // y = A*x
27    fn compute_column_means(&self) -> Vec<T>;
28    fn multiply_with_dense(&self, dense: &DMatrix<T>, result: &mut DMatrix<T>, transpose_self: bool);
29    fn multiply_with_dense_centered(&self, dense: &DMatrix<T>, result: &mut DMatrix<T>, transpose_self: bool, means: &DVector<T>);
30
31    fn multiply_transposed_by_dense(&self, q: &DMatrix<T>, result: &mut DMatrix<T>);
32    fn multiply_transposed_by_dense_centered(&self, q: &DMatrix<T>, result: &mut DMatrix<T>, means: &DVector<T>);
33}
34
35/// Singular Value Decomposition Components
36///
37/// # Fields
38/// - d:  Dimensionality (rank), the number of rows of both `ut`, `vt` and the length of `s`
39/// - ut: Transpose of left singular vectors, the vectors are the rows of `ut`
40/// - s:  Singular values (length `d`)
41/// - vt: Transpose of right singular vectors, the vectors are the rows of `vt`
42/// - diagnostics: Computational diagnostics
43#[derive(Debug, Clone, PartialEq)]
44pub struct SvdRec<T: Float> {
45    pub d: usize,
46    pub u: Array2<T>,
47    pub s: Array1<T>,
48    pub vt: Array2<T>,
49    pub diagnostics: Diagnostics<T>,
50}
51
52/// Computational Diagnostics
53///
54/// # Fields
55/// - non_zero:  Number of non-zeros in the matrix
56/// - dimensions: Number of dimensions attempted (bounded by matrix shape)
57/// - iterations: Number of iterations attempted (bounded by dimensions and matrix shape)
58/// - transposed:  True if the matrix was transposed internally
59/// - lanczos_steps: Number of Lanczos steps performed
60/// - ritz_values_stabilized: Number of ritz values
61/// - significant_values: Number of significant values discovered
62/// - singular_values: Number of singular values returned
63/// - end_interval: left, right end of interval containing unwanted eigenvalues
64/// - kappa: relative accuracy of ritz values acceptable as eigenvalues
65/// - random_seed: Random seed provided or the seed generated
66#[derive(Debug, Clone, PartialEq)]
67pub struct Diagnostics<T: Float> {
68    pub non_zero: usize,
69    pub dimensions: usize,
70    pub iterations: usize,
71    pub transposed: bool,
72    pub lanczos_steps: usize,
73    pub ritz_values_stabilized: usize,
74    pub significant_values: usize,
75    pub singular_values: usize,
76    pub end_interval: [T; 2],
77    pub kappa: T,
78    pub random_seed: u32,
79}
80
81pub trait SvdFloat: FloatOpsTS {
82    fn eps() -> Self;
83    fn eps34() -> Self;
84    fn compare(a: Self, b: Self) -> bool;
85}
86
87impl SvdFloat for f32 {
88    fn eps() -> Self {
89        f32::EPSILON
90    }
91
92    fn eps34() -> Self {
93        f32::EPSILON.powf(0.75)
94    }
95
96    fn compare(a: Self, b: Self) -> bool {
97        (b - a).abs() < f32::EPSILON
98    }
99}
100
101impl SvdFloat for f64 {
102    fn eps() -> Self {
103        f64::EPSILON
104    }
105
106    fn eps34() -> Self {
107        f64::EPSILON.powf(0.75)
108    }
109
110    fn compare(a: Self, b: Self) -> bool {
111        (b - a).abs() < f64::EPSILON
112    }
113}
single_svdlib/utils.rs

single_svdlib/
utils.rs
