scirs2_sparse/linalg/eigen/
mod.rs

1//! Eigenvalue algorithms for sparse matrices
2//!
3//! This module provides a comprehensive set of eigenvalue solvers for sparse matrices,
4//! organized by algorithm type and matrix properties.
5
6pub mod general;
7pub mod generalized;
8pub mod lanczos;
9pub mod power_iteration;
10pub mod symmetric;
11
12// Re-export main types and functions for convenience
13pub use general::eigs;
14pub use generalized::{eigsh_generalized, eigsh_generalized_enhanced};
15pub use lanczos::{lanczos, EigenResult, LanczosOptions};
16pub use power_iteration::{power_iteration, PowerIterationOptions};
17pub use symmetric::{eigsh, eigsh_shift_invert, eigsh_shift_invert_enhanced};
18
19// Common eigenvalue result type
20
21// Compatibility types for missing imports
22#[derive(Debug, Clone)]
23pub struct ArpackOptions {
24    pub max_iter: usize,
25    pub tol: f64,
26}
27
28#[derive(Debug, Clone, Copy)]
29pub enum EigenvalueMethod {
30    Lanczos,
31    PowerIteration,
32    Arpack,
33}
34
35#[derive(Debug, Clone, Copy)]
36pub enum EigenvalueMode {
37    Normal,
38    ShiftInvert,
39    Buckling,
40}
41
42/// Eigenvalue solver configuration
43#[derive(Debug, Clone)]
44#[allow(dead_code)]
45pub struct EigenSolverConfig {
46    /// Maximum number of iterations
47    pub max_iter: usize,
48    /// Convergence tolerance
49    pub tolerance: f64,
50    /// Number of eigenvalues to compute
51    pub num_eigenvalues: usize,
52    /// Whether to compute eigenvectors
53    pub compute_eigenvectors: bool,
54    /// Which eigenvalues to compute ("LA", "SA", "LM", "SM")
55    pub which: String,
56}
57
58impl Default for EigenSolverConfig {
59    fn default() -> Self {
60        Self {
61            max_iter: 1000,
62            tolerance: 1e-8,
63            num_eigenvalues: 6,
64            compute_eigenvectors: true,
65            which: "LA".to_string(),
66        }
67    }
68}
69
70impl EigenSolverConfig {
71    /// Create a new eigenvalue solver configuration
72    pub fn new() -> Self {
73        Self::default()
74    }
75
76    /// Set maximum number of iterations
77    pub fn with_max_iter(mut self, max_iter: usize) -> Self {
78        self.max_iter = max_iter;
79        self
80    }
81
82    /// Set convergence tolerance
83    pub fn with_tolerance(mut self, tolerance: f64) -> Self {
84        self.tolerance = tolerance;
85        self
86    }
87
88    /// Set number of eigenvalues to compute
89    pub fn with_num_eigenvalues(mut self, num_eigenvalues: usize) -> Self {
90        self.num_eigenvalues = num_eigenvalues;
91        self
92    }
93
94    /// Set whether to compute eigenvectors
95    pub fn with_compute_eigenvectors(mut self, compute_eigenvectors: bool) -> Self {
96        self.compute_eigenvectors = compute_eigenvectors;
97        self
98    }
99
100    /// Set which eigenvalues to compute
101    pub fn with_which(mut self, which: &str) -> Self {
102        self.which = which.to_string();
103        self
104    }
105
106    /// Convert to LanczosOptions
107    pub fn to_lanczos_options(&self) -> LanczosOptions {
108        LanczosOptions {
109            max_iter: self.max_iter,
110            max_subspace_size: (self.num_eigenvalues * 2 + 10).max(20),
111            tol: self.tolerance,
112            numeigenvalues: self.num_eigenvalues,
113            compute_eigenvectors: self.compute_eigenvectors,
114        }
115    }
116
117    /// Convert to PowerIterationOptions
118    pub fn to_power_iteration_options(&self) -> PowerIterationOptions {
119        PowerIterationOptions {
120            max_iter: self.max_iter,
121            tol: self.tolerance,
122            normalize: true,
123        }
124    }
125}
126
127/// Eigenvalue solver selection strategy
128#[derive(Debug, Clone, Copy, PartialEq, Eq)]
129pub enum EigenSolverStrategy {
130    /// Power iteration for single largest eigenvalue
131    PowerIteration,
132    /// Lanczos algorithm for multiple eigenvalues
133    Lanczos,
134    /// Symmetric eigenvalue solver (eigsh)
135    Symmetric,
136    /// Shift-invert mode for eigenvalues near a target
137    ShiftInvert,
138    /// Automatic selection based on problem characteristics
139    Auto,
140}
141
142impl Default for EigenSolverStrategy {
143    fn default() -> Self {
144        Self::Auto
145    }
146}
147
148/// Unified eigenvalue solver interface
149///
150/// This function provides a unified interface to various eigenvalue solvers,
151/// automatically selecting the best algorithm based on the problem characteristics.
152pub fn solve_eigenvalues<T>(
153    matrix: &crate::sym_csr::SymCsrMatrix<T>,
154    config: &EigenSolverConfig,
155    strategy: EigenSolverStrategy,
156) -> crate::error::SparseResult<EigenResult<T>>
157where
158    T: scirs2_core::numeric::Float
159        + scirs2_core::SparseElement
160        + std::fmt::Debug
161        + Copy
162        + std::ops::Add<Output = T>
163        + std::ops::Sub<Output = T>
164        + std::ops::Mul<Output = T>
165        + std::ops::Div<Output = T>
166        + std::iter::Sum
167        + scirs2_core::simd_ops::SimdUnifiedOps
168        + Send
169        + Sync
170        + 'static,
171{
172    let (n, _) = matrix.shape();
173
174    let selected_strategy = match strategy {
175        EigenSolverStrategy::Auto => {
176            // Automatic selection based on problem characteristics
177            if config.num_eigenvalues == 1 && config.which == "LA" {
178                EigenSolverStrategy::PowerIteration
179            } else if n < 1000 {
180                EigenSolverStrategy::Lanczos
181            } else {
182                EigenSolverStrategy::Symmetric
183            }
184        }
185        other => other,
186    };
187
188    match selected_strategy {
189        EigenSolverStrategy::PowerIteration => {
190            let opts = config.to_power_iteration_options();
191            let power_result = power_iteration::power_iteration(matrix, &opts, None)?;
192            // Convert power iteration result to common EigenResult format
193            Ok(EigenResult {
194                eigenvalues: power_result.eigenvalues,
195                eigenvectors: power_result.eigenvectors,
196                iterations: power_result.iterations,
197                residuals: power_result.residuals,
198                converged: power_result.converged,
199            })
200        }
201        EigenSolverStrategy::Lanczos => {
202            let opts = config.to_lanczos_options();
203            let result = lanczos::lanczos(matrix, &opts, None)?;
204            Ok(result)
205        }
206        EigenSolverStrategy::Symmetric => {
207            let opts = config.to_lanczos_options();
208            let result = symmetric::eigsh(
209                matrix,
210                Some(config.num_eigenvalues),
211                Some(&config.which),
212                Some(opts),
213            )?;
214            Ok(result)
215        }
216        EigenSolverStrategy::ShiftInvert => {
217            // For shift-invert, we need a target value - use 0.0 as default
218            let opts = config.to_lanczos_options();
219            let result = symmetric::eigsh_shift_invert(
220                matrix,
221                T::sparse_zero(),
222                Some(config.num_eigenvalues),
223                Some(&config.which),
224                Some(opts),
225            )?;
226            Ok(result)
227        }
228        EigenSolverStrategy::Auto => {
229            // This case is handled above, but we include it for completeness
230            unreachable!("Auto strategy should have been resolved")
231        }
232    }
233}
234
235#[cfg(test)]
236mod tests {
237    use super::*;
238    use crate::sym_csr::SymCsrMatrix;
239
240    #[test]
241    fn test_eigen_solver_config() {
242        let config = EigenSolverConfig::new()
243            .with_max_iter(500)
244            .with_tolerance(1e-10)
245            .with_num_eigenvalues(3)
246            .with_which("SA");
247
248        assert_eq!(config.max_iter, 500);
249        assert_eq!(config.tolerance, 1e-10);
250        assert_eq!(config.num_eigenvalues, 3);
251        assert_eq!(config.which, "SA");
252    }
253
254    #[test]
255    fn test_unified_solver_auto() {
256        // Create a simple 2x2 symmetric matrix
257        // Matrix: [[2, 1], [1, 2]] stored as lower: [[2], [1, 2]]
258        let data = vec![2.0, 1.0, 2.0];
259        let indptr = vec![0, 1, 3];
260        let indices = vec![0, 0, 1];
261        let matrix = SymCsrMatrix::new(data, indptr, indices, (2, 2)).unwrap();
262
263        let config = EigenSolverConfig::new().with_num_eigenvalues(1);
264        let result = solve_eigenvalues(&matrix, &config, EigenSolverStrategy::Auto).unwrap();
265
266        assert!(result.converged);
267        assert_eq!(result.eigenvalues.len(), 1);
268    }
269
270    #[test]
271    fn test_unified_solver_power_iteration() {
272        // Matrix: [[3, 1], [1, 2]] stored as lower: [[3], [1, 2]]
273        let data = vec![3.0, 1.0, 2.0];
274        let indptr = vec![0, 1, 3];
275        let indices = vec![0, 0, 1];
276        let matrix = SymCsrMatrix::new(data, indptr, indices, (2, 2)).unwrap();
277
278        let config = EigenSolverConfig::new().with_num_eigenvalues(1);
279        let result =
280            solve_eigenvalues(&matrix, &config, EigenSolverStrategy::PowerIteration).unwrap();
281
282        assert!(result.converged);
283        assert_eq!(result.eigenvalues.len(), 1);
284    }
285
286    #[test]
287    fn test_unified_solver_lanczos() {
288        // Matrix: [[4, 2], [2, 3]] stored as lower: [[4], [2, 3]]
289        let data = vec![4.0, 2.0, 3.0];
290        let indptr = vec![0, 1, 3];
291        let indices = vec![0, 0, 1];
292        let matrix = SymCsrMatrix::new(data, indptr, indices, (2, 2)).unwrap();
293
294        let config = EigenSolverConfig::new().with_num_eigenvalues(2);
295        let result = solve_eigenvalues(&matrix, &config, EigenSolverStrategy::Lanczos).unwrap();
296
297        assert!(result.converged);
298        assert!(!result.eigenvalues.is_empty());
299    }
300}
scirs2_sparse/linalg/eigen/mod.rs

scirs2_sparse/linalg/eigen/
mod.rs
