Reorganized structure into multiple modules

Author: Ian

src/masked.rs src/laczos/masked.rssrc/masked.rs renamed to src/laczos/masked.rs

@@ -76,7 +76,7 @@ impl<'a, T: Float + AddAssign + Sync + Send> SMat<T> for MaskedCSRMatrix<'a, T>
 
         for i in 0..self.matrix.nrows() {
             for j in major_offsets[i]..major_offsets[i + 1] {
-                let col = minor_indices[j]; 
+                let col = minor_indices[j];
                 if self.column_mask[col] {
                     count += 1;
                 }
@@ -213,7 +213,7 @@ impl<'a, T: Float + AddAssign + Sync + Send> SMat<T> for MaskedCSRMatrix<'a, T>
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::{svd, SMat};
+    use crate::{SMat};
     use nalgebra_sparse::{coo::CooMatrix, csr::CsrMatrix};
     use rand::rngs::StdRng;
     use rand::{Rng, SeedableRng};
@@ -243,7 +243,7 @@ mod tests {
         assert_eq!(masked.nnz(), 6); // Only entries in the selected columns
 
         // Test SVD on the masked matrix
-        let svd_result = svd(&masked);
+        let svd_result = crate::laczos::svd(&masked);
         assert!(svd_result.is_ok());
     }
 
@@ -304,8 +304,8 @@ mod tests {
         assert_eq!(masked_matrix.nnz(), physical_csr.nnz());
 
         // Perform SVD on both
-        let svd_masked = svd(&masked_matrix).unwrap();
-        let svd_physical = svd(&physical_csr).unwrap();
+        let svd_masked = crate::laczos::svd(&masked_matrix).unwrap();
+        let svd_physical = crate::laczos::svd(&physical_csr).unwrap();
 
         // Compare SVD results - they should be very close but not exactly the same
         // due to potential differences in numerical computation

src/new.rs src/laczos/mod.rssrc/new.rs renamed to src/laczos/mod.rs

@@ -1,5 +1,8 @@
+pub mod masked;
+
 use crate::error::SvdLibError;
-use ndarray::{Array, Array1, Array2};
+use crate::{Diagnostics, SMat, SvdFloat, SvdRec};
+use ndarray::{Array, Array2};
 use num_traits::real::Real;
 use num_traits::{Float, FromPrimitive, One, Zero};
 use rand::rngs::StdRng;
@@ -12,106 +15,7 @@ use std::iter::Sum;
 use std::mem;
 use std::ops::{AddAssign, MulAssign, Neg, SubAssign};
 
-pub trait SMat<T: Float>: Sync {
-    fn nrows(&self) -> usize;
-    fn ncols(&self) -> usize;
-    fn nnz(&self) -> usize;
-    fn svd_opa(&self, x: &[T], y: &mut [T], transposed: bool); // y = A*x
-}
-
-/// Singular Value Decomposition Components
-///
-/// # Fields
-/// - d:  Dimensionality (rank), the number of rows of both `ut`, `vt` and the length of `s`
-/// - ut: Transpose of left singular vectors, the vectors are the rows of `ut`
-/// - s:  Singular values (length `d`)
-/// - vt: Transpose of right singular vectors, the vectors are the rows of `vt`
-/// - diagnostics: Computational diagnostics
-#[derive(Debug, Clone, PartialEq)]
-pub struct SvdRec<T: Float> {
-    pub d: usize,
-    pub ut: Array2<T>,
-    pub s: Array1<T>,
-    pub vt: Array2<T>,
-    pub diagnostics: Diagnostics<T>,
-}
-
-/// Computational Diagnostics
-///
-/// # Fields
-/// - non_zero:  Number of non-zeros in the matrix
-/// - dimensions: Number of dimensions attempted (bounded by matrix shape)
-/// - iterations: Number of iterations attempted (bounded by dimensions and matrix shape)
-/// - transposed:  True if the matrix was transposed internally
-/// - lanczos_steps: Number of Lanczos steps performed
-/// - ritz_values_stabilized: Number of ritz values
-/// - significant_values: Number of significant values discovered
-/// - singular_values: Number of singular values returned
-/// - end_interval: left, right end of interval containing unwanted eigenvalues
-/// - kappa: relative accuracy of ritz values acceptable as eigenvalues
-/// - random_seed: Random seed provided or the seed generated
-#[derive(Debug, Clone, PartialEq)]
-pub struct Diagnostics<T: Float> {
-    pub non_zero: usize,
-    pub dimensions: usize,
-    pub iterations: usize,
-    pub transposed: bool,
-    pub lanczos_steps: usize,
-    pub ritz_values_stabilized: usize,
-    pub significant_values: usize,
-    pub singular_values: usize,
-    pub end_interval: [T; 2],
-    pub kappa: T,
-    pub random_seed: u32,
-}
-
 /// Trait for floating point types that can be used with the SVD algorithm
-pub trait SvdFloat:
-    Float
-    + FromPrimitive
-    + Debug
-    + Send
-    + Sync
-    + Zero
-    + One
-    + AddAssign
-    + SubAssign
-    + MulAssign
-    + Neg<Output = Self>
-    + Sum
-{
-    fn eps() -> Self;
-    fn eps34() -> Self;
-    fn compare(a: Self, b: Self) -> bool;
-}
-
-impl SvdFloat for f32 {
-    fn eps() -> Self {
-        f32::EPSILON
-    }
-
-    fn eps34() -> Self {
-        f32::EPSILON.powf(0.75)
-    }
-
-    fn compare(a: Self, b: Self) -> bool {
-        (b - a).abs() < f32::EPSILON
-    }
-}
-
-impl SvdFloat for f64 {
-    fn eps() -> Self {
-        f64::EPSILON
-    }
-
-    fn eps34() -> Self {
-        f64::EPSILON.powf(0.75)
-    }
-
-    fn compare(a: Self, b: Self) -> bool {
-        (b - a).abs() < f64::EPSILON
-    }
-}
 
 /// SVD at full dimensionality, calls `svdLAS2` with the highlighted defaults
 ///
@@ -1229,7 +1133,6 @@ fn ritvec<T: SvdFloat>(
     let significant_indices: Vec<usize> = (0..js)
         .into_par_iter()
         .filter(|&k| {
-            // Adaptive error bound check using relative tolerance
             let relative_bound =
                 adaptive_kappa * wrk.ritz[k].abs().max(max_eigenvalue * adaptive_eps);
             wrk.bnd[k] <= relative_bound && k + 1 > js - neig
@@ -1242,19 +1145,17 @@ fn ritvec<T: SvdFloat>(
         .into_par_iter()
         .map(|k| {
             let mut vec = vec![T::zero(); wrk.ncols];
-            let mut idx = (jsq - js) + k + 1;
 
             for i in 0..js {
-                idx -= js;
-                // Non-zero check with adaptive threshold
+                let idx = k * js + i;
+
                 if s[idx].abs() > adaptive_eps {
                     for (j, item) in store_vectors[i].iter().enumerate().take(wrk.ncols) {
                         vec[j] += s[idx] * *item;
                     }
                 }
             }
 
-            // Return with position index (for proper ordering)
             (k, vec)
         })
         .collect();

src/lib.rs

@@ -1,14 +1,13 @@
 pub mod legacy;
 pub mod error;
-mod new;
-mod masked;
 pub(crate) mod utils;
 
-pub(crate) mod randomized;
+pub mod randomized;
+
+pub mod laczos;
+
+pub use utils::*;
 
-pub use new::*;
-pub use masked::*;
-pub use randomized::{randomized_svd, PowerIterationNormalizer};
 
 #[cfg(test)]
 mod simple_comparison_tests {
@@ -19,7 +18,6 @@ mod simple_comparison_tests {
     use rand::{Rng, SeedableRng};
     use rand::rngs::StdRng;
     use rayon::ThreadPoolBuilder;
-    use crate::randomized::randomized_svd;
 
     fn create_sparse_matrix(rows: usize, cols: usize, density: f64) -> nalgebra_sparse::coo::CooMatrix<f64> {
         use rand::{rngs::StdRng, Rng, SeedableRng};
@@ -77,7 +75,7 @@ mod simple_comparison_tests {
 
         // Run both implementations with the same seed for deterministic behavior
         let seed = 42;
-        let current_result = svd_dim_seed(&test_matrix, 0, seed).unwrap();
+        let current_result = laczos::svd_dim_seed(&test_matrix, 0, seed).unwrap();
         let legacy_result = legacy::svd_dim_seed(&test_matrix, 0, seed).unwrap();
 
         // Compare dimensions
@@ -129,12 +127,12 @@ mod simple_comparison_tests {
         let csr = CsrMatrix::from(&coo);
 
         // Calculate SVD using original method
-        let legacy_svd = svd_dim_seed(&csr, 0, seed as u32).unwrap();
+        let legacy_svd = laczos::svd_dim_seed(&csr, 0, seed as u32).unwrap();
 
         // Calculate SVD using our masked method (using all columns)
         let mask = vec![true; ncols];
-        let masked_matrix = MaskedCSRMatrix::new(&csr, mask);
-        let current_svd = svd_dim_seed(&masked_matrix, 0, seed as u32).unwrap();
+        let masked_matrix = laczos::masked::MaskedCSRMatrix::new(&csr, mask);
+        let current_svd = laczos::svd_dim_seed(&masked_matrix, 0, seed as u32).unwrap();
 
         // Compare with relative tolerance
         let rel_tol = 1e-3;  // 0.1% relative tolerance
@@ -161,13 +159,12 @@ mod simple_comparison_tests {
         let test_matrix = create_sparse_matrix(100, 100, 0.0098); // 0.98% non-zeros
 
         // Should no longer fail with convergence error
-        let result = svd_dim_seed(&test_matrix, 50, 42);
+        let result = laczos::svd_dim_seed(&test_matrix, 50, 42);
         assert!(result.is_ok(), "{}", format!("SVD failed on 99.02% sparse matrix, {:?}", result.err().unwrap()));
     }
 
     #[test]
     fn test_random_svd_computation() {
-        use crate::{randomized_svd, PowerIterationNormalizer};
 
         // Create a matrix with high sparsity (99%)
         let test_matrix = create_sparse_matrix(1000, 250, 0.01); // 1% non-zeros
@@ -176,12 +173,12 @@ mod simple_comparison_tests {
         let csr = CsrMatrix::from(&test_matrix);
 
         // Run randomized SVD with reasonable defaults for a sparse matrix
-        let result = randomized_svd(
+        let result = randomized::randomized_svd(
             &csr,
             50,                              // target rank
             10,                              // oversampling parameter
             3,                               // power iterations
-            PowerIterationNormalizer::QR,    // use QR normalization
+            randomized::PowerIterationNormalizer::QR,    // use QR normalization
             Some(42),                        // random seed
         );
 
@@ -219,7 +216,6 @@ mod simple_comparison_tests {
 
     #[test]
     fn test_randomized_svd_very_large_sparse_matrix() {
-        use crate::{randomized_svd, PowerIterationNormalizer};
 
         // Create a very large matrix with high sparsity (99%)
         let test_matrix = create_sparse_matrix(100000, 2500, 0.01); // 1% non-zeros
@@ -230,12 +226,12 @@ mod simple_comparison_tests {
         // Run randomized SVD with reasonable defaults for a sparse matrix
         let threadpool = ThreadPoolBuilder::new().num_threads(10).build().unwrap();
         let result = threadpool.install(|| {
-             randomized_svd(
+            randomized::randomized_svd(
                 &csr,
                 50,                              // target rank
                 10,                              // oversampling parameter
-                2,                               // power iterations
-                PowerIterationNormalizer::QR,    // use QR normalization
+                7,                               // power iterations
+                randomized::PowerIterationNormalizer::QR,    // use QR normalization
                 Some(42),                        // random seed
             )
         });
@@ -251,7 +247,6 @@ mod simple_comparison_tests {
 
     #[test]
     fn test_randomized_svd_small_sparse_matrix() {
-        use crate::{randomized_svd, PowerIterationNormalizer};
 
         // Create a very large matrix with high sparsity (99%)
         let test_matrix = create_sparse_matrix(1000, 250, 0.01); // 1% non-zeros
@@ -262,12 +257,12 @@ mod simple_comparison_tests {
         // Run randomized SVD with reasonable defaults for a sparse matrix
         let threadpool = ThreadPoolBuilder::new().num_threads(10).build().unwrap();
         let result = threadpool.install(|| {
-            randomized_svd(
+            randomized::randomized_svd(
                 &csr,
                 50,                              // target rank
                 10,                              // oversampling parameter
                 2,                               // power iterations
-                PowerIterationNormalizer::QR,    // use QR normalization
+                randomized::PowerIterationNormalizer::QR,    // use QR normalization
                 Some(42),                        // random seed
             )
         });