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

scirs2_sparse/linalg/eigen/
mod.rs
