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        + std::fmt::Debug
160        + Copy
161        + std::ops::Add<Output = T>
162        + std::ops::Sub<Output = T>
163        + std::ops::Mul<Output = T>
164        + std::ops::Div<Output = T>
165        + std::iter::Sum
166        + scirs2_core::simd_ops::SimdUnifiedOps
167        + Send
168        + Sync
169        + 'static,
170{
171    let (n, _) = matrix.shape();
172
173    let selected_strategy = match strategy {
174        EigenSolverStrategy::Auto => {
175            // Automatic selection based on problem characteristics
176            if config.num_eigenvalues == 1 && config.which == "LA" {
177                EigenSolverStrategy::PowerIteration
178            } else if n < 1000 {
179                EigenSolverStrategy::Lanczos
180            } else {
181                EigenSolverStrategy::Symmetric
182            }
183        }
184        other => other,
185    };
186
187    match selected_strategy {
188        EigenSolverStrategy::PowerIteration => {
189            let opts = config.to_power_iteration_options();
190            let power_result = power_iteration::power_iteration(matrix, &opts, None)?;
191            // Convert power iteration result to common EigenResult format
192            Ok(EigenResult {
193                eigenvalues: power_result.eigenvalues,
194                eigenvectors: power_result.eigenvectors,
195                iterations: power_result.iterations,
196                residuals: power_result.residuals,
197                converged: power_result.converged,
198            })
199        }
200        EigenSolverStrategy::Lanczos => {
201            let opts = config.to_lanczos_options();
202            let result = lanczos::lanczos(matrix, &opts, None)?;
203            Ok(result)
204        }
205        EigenSolverStrategy::Symmetric => {
206            let opts = config.to_lanczos_options();
207            let result = symmetric::eigsh(
208                matrix,
209                Some(config.num_eigenvalues),
210                Some(&config.which),
211                Some(opts),
212            )?;
213            Ok(result)
214        }
215        EigenSolverStrategy::ShiftInvert => {
216            // For shift-invert, we need a target value - use 0.0 as default
217            let opts = config.to_lanczos_options();
218            let result = symmetric::eigsh_shift_invert(
219                matrix,
220                T::zero(),
221                Some(config.num_eigenvalues),
222                Some(&config.which),
223                Some(opts),
224            )?;
225            Ok(result)
226        }
227        EigenSolverStrategy::Auto => {
228            // This case is handled above, but we include it for completeness
229            unreachable!("Auto strategy should have been resolved")
230        }
231    }
232}
233
234#[cfg(test)]
235mod tests {
236    use super::*;
237    use crate::sym_csr::SymCsrMatrix;
238
239    #[test]
240    fn test_eigen_solver_config() {
241        let config = EigenSolverConfig::new()
242            .with_max_iter(500)
243            .with_tolerance(1e-10)
244            .with_num_eigenvalues(3)
245            .with_which("SA");
246
247        assert_eq!(config.max_iter, 500);
248        assert_eq!(config.tolerance, 1e-10);
249        assert_eq!(config.num_eigenvalues, 3);
250        assert_eq!(config.which, "SA");
251    }
252
253    #[test]
254    fn test_unified_solver_auto() {
255        // Create a simple 2x2 symmetric matrix
256        // Matrix: [[2, 1], [1, 2]] stored as lower: [[2], [1, 2]]
257        let data = vec![2.0, 1.0, 2.0];
258        let indptr = vec![0, 1, 3];
259        let indices = vec![0, 0, 1];
260        let matrix = SymCsrMatrix::new(data, indptr, indices, (2, 2)).unwrap();
261
262        let config = EigenSolverConfig::new().with_num_eigenvalues(1);
263        let result = solve_eigenvalues(&matrix, &config, EigenSolverStrategy::Auto).unwrap();
264
265        assert!(result.converged);
266        assert_eq!(result.eigenvalues.len(), 1);
267    }
268
269    #[test]
270    fn test_unified_solver_power_iteration() {
271        // Matrix: [[3, 1], [1, 2]] stored as lower: [[3], [1, 2]]
272        let data = vec![3.0, 1.0, 2.0];
273        let indptr = vec![0, 1, 3];
274        let indices = vec![0, 0, 1];
275        let matrix = SymCsrMatrix::new(data, indptr, indices, (2, 2)).unwrap();
276
277        let config = EigenSolverConfig::new().with_num_eigenvalues(1);
278        let result =
279            solve_eigenvalues(&matrix, &config, EigenSolverStrategy::PowerIteration).unwrap();
280
281        assert!(result.converged);
282        assert_eq!(result.eigenvalues.len(), 1);
283    }
284
285    #[test]
286    fn test_unified_solver_lanczos() {
287        // Matrix: [[4, 2], [2, 3]] stored as lower: [[4], [2, 3]]
288        let data = vec![4.0, 2.0, 3.0];
289        let indptr = vec![0, 1, 3];
290        let indices = vec![0, 0, 1];
291        let matrix = SymCsrMatrix::new(data, indptr, indices, (2, 2)).unwrap();
292
293        let config = EigenSolverConfig::new().with_num_eigenvalues(2);
294        let result = solve_eigenvalues(&matrix, &config, EigenSolverStrategy::Lanczos).unwrap();
295
296        assert!(result.converged);
297        assert!(!result.eigenvalues.is_empty());
298    }
299}
scirs2_sparse/linalg/eigen/mod.rs

scirs2_sparse/linalg/eigen/
mod.rs
