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