scirs2_sparse/
lib.rs

1#![allow(clippy::manual_div_ceil)]
2#![allow(clippy::needless_return)]
3#![allow(clippy::manual_ok_err)]
4#![allow(clippy::needless_range_loop)]
5#![allow(clippy::while_let_loop)]
6#![allow(clippy::vec_init_then_push)]
7#![allow(clippy::should_implement_trait)]
8#![allow(clippy::only_used_in_recursion)]
9#![allow(clippy::manual_slice_fill)]
10#![allow(dead_code)]
11//! # SciRS2 Sparse - Sparse Matrix Operations
12//!
13//! **scirs2-sparse** provides comprehensive sparse matrix formats and operations modeled after SciPy's
14//! `sparse` module, offering CSR, CSC, COO, DOK, LIL, DIA, BSR formats with efficient algorithms
15//! for large-scale sparse linear algebra, eigenvalue problems, and graph operations.
16//!
17//! ## 🎯 Key Features
18//!
19//! - **SciPy Compatibility**: Drop-in replacement for `scipy.sparse` classes
20//! - **Multiple Formats**: CSR, CSC, COO, DOK, LIL, DIA, BSR with easy conversion
21//! - **Efficient Operations**: Sparse matrix-vector/matrix multiplication
22//! - **Linear Solvers**: Direct (LU, Cholesky) and iterative (CG, GMRES) solvers
23//! - **Eigenvalue Solvers**: ARPACK-based sparse eigenvalue computation
24//! - **Array API**: Modern NumPy-compatible array interface (recommended)
25//!
26//! ## 📦 Module Overview
27//!
28//! | SciRS2 Format | SciPy Equivalent | Description |
29//! |---------------|------------------|-------------|
30//! | `CsrArray` | `scipy.sparse.csr_array` | Compressed Sparse Row (efficient row slicing) |
31//! | `CscArray` | `scipy.sparse.csc_array` | Compressed Sparse Column (efficient column slicing) |
32//! | `CooArray` | `scipy.sparse.coo_array` | Coordinate format (efficient construction) |
33//! | `DokArray` | `scipy.sparse.dok_array` | Dictionary of Keys (efficient element access) |
34//! | `LilArray` | `scipy.sparse.lil_array` | List of Lists (efficient incremental construction) |
35//! | `DiaArray` | `scipy.sparse.dia_array` | Diagonal format (efficient banded matrices) |
36//! | `BsrArray` | `scipy.sparse.bsr_array` | Block Sparse Row (efficient block operations) |
37//!
38//! ## 🚀 Quick Start
39//!
40//! ```toml
41//! [dependencies]
42//! scirs2-sparse = "0.4.0"
43//! ```
44//!
45//! ```rust
46//! use scirs2_sparse::csr_array::CsrArray;
47//!
48//! // Create sparse matrix from triplets (row, col, value)
49//! let rows = vec![0, 0, 1, 2, 2];
50//! let cols = vec![0, 2, 2, 0, 1];
51//! let data = vec![1.0, 2.0, 3.0, 4.0, 5.0];
52//! let sparse = CsrArray::from_triplets(&rows, &cols, &data, (3, 3), false).expect("Operation failed");
53//! ```
54//!
55//! ## 🔒 Version: 0.2.0 (February 8, 2026)
56//!
57//! ## Matrix vs. Array API
58//!
59//! This module provides both a matrix-based API and an array-based API,
60//! following SciPy's transition to a more NumPy-compatible array interface.
61//!
62//! When using the array interface (e.g., `CsrArray`), please note that:
63//!
64//! - `*` performs element-wise multiplication, not matrix multiplication
65//! - Use `dot()` method for matrix multiplication
66//! - Operations like `sum` produce arrays, not matrices
67//! - Array-style slicing operations return scalars, 1D, or 2D arrays
68//!
69//! For new code, we recommend using the array interface, which is more consistent
70//! with the rest of the numerical ecosystem.
71//!
72//! ## Examples
73//!
74//! ### Matrix API (Legacy)
75//!
76//! ```
77//! use scirs2_sparse::csr::CsrMatrix;
78//!
79//! // Create a sparse matrix in CSR format
80//! let rows = vec![0, 0, 1, 2, 2];
81//! let cols = vec![0, 2, 2, 0, 1];
82//! let data = vec![1.0, 2.0, 3.0, 4.0, 5.0];
83//! let shape = (3, 3);
84//!
85//! let matrix = CsrMatrix::new(data, rows, cols, shape).expect("Operation failed");
86//! ```
87//!
88//! ### Array API (Recommended)
89//!
90//! ```
91//! use scirs2_sparse::csr_array::CsrArray;
92//!
93//! // Create a sparse array in CSR format
94//! let rows = vec![0, 0, 1, 2, 2];
95//! let cols = vec![0, 2, 2, 0, 1];
96//! let data = vec![1.0, 2.0, 3.0, 4.0, 5.0];
97//! let shape = (3, 3);
98//!
99//! // From triplets (COO-like construction)
100//! let array = CsrArray::from_triplets(&rows, &cols, &data, shape, false).expect("Operation failed");
101//!
102//! // Or directly from CSR components
103//! // let array = CsrArray::new(...);
104//! ```
105
106// Export error types
107pub mod error;
108pub use error::{SparseError, SparseResult};
109
110// Base trait for sparse arrays
111pub mod sparray;
112pub use sparray::{is_sparse, SparseArray, SparseSum};
113
114// Trait for symmetric sparse arrays
115pub mod sym_sparray;
116pub use sym_sparray::SymSparseArray;
117
118// No spatial module in sparse
119
120// Array API (recommended)
121pub mod csr_array;
122pub use csr_array::CsrArray;
123
124pub mod csc_array;
125pub use csc_array::CscArray;
126
127pub mod coo_array;
128pub use coo_array::CooArray;
129
130pub mod dok_array;
131pub use dok_array::DokArray;
132
133pub mod lil_array;
134pub use lil_array::LilArray;
135
136pub mod dia_array;
137pub use dia_array::DiaArray;
138
139pub mod bsr_array;
140pub use bsr_array::BsrArray;
141
142pub mod banded_array;
143pub use banded_array::BandedArray;
144
145// Symmetric array formats
146pub mod sym_csr;
147pub use sym_csr::{SymCsrArray, SymCsrMatrix};
148
149pub mod sym_coo;
150pub use sym_coo::{SymCooArray, SymCooMatrix};
151
152// Legacy matrix formats
153pub mod csr;
154pub use csr::CsrMatrix;
155
156pub mod csc;
157pub use csc::CscMatrix;
158
159pub mod coo;
160pub use coo::CooMatrix;
161
162pub mod dok;
163pub use dok::DokMatrix;
164
165pub mod lil;
166pub use lil::LilMatrix;
167
168pub mod dia;
169pub use dia::DiaMatrix;
170
171pub mod bsr;
172pub use bsr::BsrMatrix;
173
174pub mod banded;
175pub use banded::BandedMatrix;
176
177// Utility functions
178pub mod utils;
179
180// Linear algebra with sparse matrices
181pub mod linalg;
182// Re-export the main functions from the reorganized linalg module
183pub use linalg::{
184    // Functions from solvers
185    add,
186    // Functions from iterative
187    bicg,
188    bicgstab,
189    cg,
190    cholesky_decomposition,
191    // Enhanced operators
192    convolution_operator,
193    diag_matrix,
194    eigs,
195    eigsh,
196    enhanced_add,
197    enhanced_diagonal,
198    enhanced_scale,
199    enhanced_subtract,
200    expm,
201    // Functions from matfuncs
202    expm_multiply,
203    eye,
204    finite_difference_operator,
205    // GCROT solver
206    gcrot,
207    gmres,
208    incomplete_cholesky,
209    incomplete_lu,
210    inv,
211    // IRAM eigenvalue solver (v0.3.0)
212    iram,
213    iram_shift_invert,
214    lanczos,
215    // Decomposition functions
216    lu_decomposition,
217    matmul,
218    matrix_power,
219    multiply,
220    norm,
221    onenormest,
222    // Eigenvalue functions
223    power_iteration,
224    qr_decomposition,
225    // Specialized solvers (v0.2.0)
226    solve_arrow_matrix,
227    solve_banded_system,
228    solve_block_2x2,
229    solve_kronecker_system,
230    solve_saddle_point,
231    sparse_direct_solve,
232    sparse_lstsq,
233    spsolve,
234    svd_truncated,
235    // SVD functions
236    svds,
237    // TFQMR solver
238    tfqmr,
239    // IRAM (Arnoldi) eigenvalue solver (v0.3.0)
240    ArnoldiConfig,
241    ArpackOptions,
242    // Interfaces
243    AsLinearOperator,
244    // Types from iterative
245    BiCGOptions,
246    BiCGSTABOptions,
247    BiCGSTABResult,
248    // Enhanced operator types
249    BoundaryCondition,
250    CGOptions,
251    CGSOptions,
252    CGSResult,
253    CholeskyResult,
254    ConvolutionMode,
255    ConvolutionOperator,
256    // Operator types
257    DiagonalOperator,
258    EigenResult,
259    EigenvalueMethod,
260    EnhancedDiagonalOperator,
261    EnhancedDifferenceOperator,
262    EnhancedOperatorOptions,
263    EnhancedScaledOperator,
264    EnhancedSumOperator,
265    FiniteDifferenceOperator,
266    GCROTOptions,
267    GCROTResult,
268    GMRESOptions,
269    ICOptions,
270    // Preconditioners
271    ILU0Preconditioner,
272    ILUOptions,
273    IdentityOperator,
274    IterationResult,
275    JacobiPreconditioner,
276    // Decomposition types
277    LUResult,
278    LanczosOptions,
279    LinearOperator,
280    // Eigenvalue types
281    PowerIterationOptions,
282    QRResult,
283    SSORPreconditioner,
284    // SVD types
285    SVDOptions,
286    SVDResult,
287    ScaledIdentityOperator,
288    TFQMROptions,
289    TFQMRResult,
290};
291
292// Format conversions
293pub mod convert;
294
295// Construction utilities
296pub mod construct;
297pub mod construct_sym;
298
299// Combining arrays
300pub mod combine;
301pub use combine::{block_diag, bmat, hstack, kron, kronsum, tril, triu, vstack};
302
303// Index dtype handling utilities
304pub mod index_dtype;
305pub use index_dtype::{can_cast_safely, get_index_dtype, safely_cast_index_arrays};
306
307// Optimized operations for symmetric sparse formats
308pub mod sym_ops;
309
310// Tensor-based sparse operations (v0.2.0)
311pub mod tensor_sparse;
312
313// Enhanced BSR format with flat storage and Block LU (v0.3.0)
314pub mod bsr_enhanced;
315pub use bsr_enhanced::{block_lu, block_lu_solve, BlockLUResult, EnhancedBsr};
316
317// Enhanced DIA format with banded operations (v0.3.0)
318pub mod dia_enhanced;
319pub use dia_enhanced::{banded_solve, dia_tridiagonal_solve, tridiagonal_solve, EnhancedDia};
320
321// Compressed Sparse Fiber (CSF) tensor format (v0.3.0)
322pub mod csf_tensor;
323pub use csf_tensor::CsfTensor;
324
325// Advanced sparse formats (v0.4.0): SELL-C-sigma, CSR5, CSF (standalone), Polynomial Preconditioners
326pub mod formats;
327pub use formats::csf::CsfTensor as CsfTensorStandalone;
328pub use formats::csr5::Csr5Matrix;
329pub use formats::poly_precond::{
330    ChebyshevPreconditioner, NeumannPreconditioner as NeumannPreconditionerPoly, PolyPrecondConfig,
331};
332pub use formats::sell::SellMatrix;
333
334// Sparse matrix utility functions (v0.3.0)
335pub mod sparse_functions;
336pub use sparse_functions::{
337    sparse_block_diag, sparse_diag_matrix, sparse_diags, sparse_eye, sparse_eye_rect,
338    sparse_hstack, sparse_kron, sparse_kronsum, sparse_random, sparse_vstack,
339};
340pub use sym_ops::{
341    sym_coo_matvec, sym_csr_matvec, sym_csr_quadratic_form, sym_csr_rank1_update, sym_csr_trace,
342};
343
344// Tensor operations (v0.2.0)
345pub use tensor_sparse::{khatri_rao_product, CPDecomposition, SparseTensor, TuckerDecomposition};
346
347// GPU-accelerated operations
348pub mod gpu;
349pub mod gpu_kernel_execution;
350pub mod gpu_ops;
351pub mod gpu_spmv_implementation;
352pub use gpu_kernel_execution::{
353    calculate_adaptive_workgroup_size, execute_spmv_kernel, execute_symmetric_spmv_kernel,
354    execute_triangular_solve_kernel, GpuKernelConfig, GpuMemoryManager as GpuKernelMemoryManager,
355    GpuPerformanceProfiler, MemoryStrategy,
356};
357pub use gpu_ops::{
358    gpu_sparse_matvec, gpu_sym_sparse_matvec, AdvancedGpuOps, GpuKernelScheduler, GpuMemoryManager,
359    GpuOptions, GpuProfiler, OptimizedGpuOps,
360};
361pub use gpu_spmv_implementation::GpuSpMV;
362
363// Memory-efficient algorithms and patterns
364pub mod memory_efficient;
365pub use memory_efficient::{
366    streaming_sparse_matvec, CacheAwareOps, MemoryPool, MemoryTracker, OutOfCoreProcessor,
367};
368
369// SIMD-accelerated operations
370pub mod simd_ops;
371pub use simd_ops::{
372    simd_csr_matvec, simd_sparse_elementwise, simd_sparse_linear_combination, simd_sparse_matmul,
373    simd_sparse_norm, simd_sparse_scale, simd_sparse_transpose, ElementwiseOp, SimdOptions,
374};
375
376// Parallel vector operations for iterative solvers
377pub mod parallel_vector_ops;
378pub use parallel_vector_ops::{
379    advanced_sparse_matvec_csr, parallel_axpy, parallel_dot, parallel_linear_combination,
380    parallel_norm2, parallel_sparse_matvec_csr, parallel_vector_add, parallel_vector_copy,
381    parallel_vector_scale, parallel_vector_sub, ParallelVectorOptions,
382};
383
384// Enhanced iterative solvers with preconditioners and sparse utilities (v0.3.0)
385pub mod iterative_solvers;
386pub use iterative_solvers::{
387    // Solvers
388    bicgstab as enhanced_bicgstab,
389    cg as enhanced_cg,
390    chebyshev,
391    // Sparse utility functions
392    estimate_spectral_radius,
393    gmres as enhanced_gmres,
394    sparse_diagonal,
395    sparse_norm,
396    sparse_trace,
397    // Preconditioners (Array1-based interface)
398    ILU0Preconditioner as EnhancedILU0Preconditioner,
399    // Configuration and result types
400    IterativeSolverConfig,
401    JacobiPreconditioner as EnhancedJacobiPreconditioner,
402    NormType,
403    Preconditioner,
404    SSORPreconditioner as EnhancedSSORPreconditioner,
405    SolverResult,
406};
407
408// LOBPCG eigensolver (v0.3.0)
409pub mod lobpcg;
410pub use lobpcg::{
411    lobpcg as lobpcg_eigensolver, lobpcg_generalised, EigenTarget, LobpcgConfig, LobpcgResult,
412};
413
414// Advanced Krylov subspace eigensolvers (v0.3.0)
415pub mod krylov;
416pub use krylov::{
417    iram as krylov_iram, thick_restart_lanczos, IramConfig, KrylovEigenResult,
418    ThickRestartLanczosConfig, WhichEigenvalues,
419};
420
421// Sparse matrix utilities (v0.3.0)
422pub mod sparse_utils;
423pub use sparse_utils::{
424    condest_1norm, permute_matrix, reverse_cuthill_mckee, sparse_add, sparse_extract_diagonal,
425    sparse_identity, sparse_kronecker, sparse_matrix_norm, sparse_matrix_trace, sparse_scale,
426    sparse_sub, sparse_transpose, spgemm, RcmResult, SparseNorm,
427};
428
429// Incomplete factorization preconditioners (v0.3.0)
430pub mod incomplete_factorizations;
431pub use incomplete_factorizations::{Ic0, Ilu0 as Ilu0Enhanced, IluK, Ilut, IlutConfig, Milu};
432
433// Sparse direct solvers (v0.3.0)
434pub mod direct_solver;
435pub use direct_solver::{
436    amd_ordering, elimination_tree, inverse_perm, nested_dissection_ordering,
437    sparse_cholesky_solve, sparse_lu_solve, symbolic_cholesky, SparseCholResult,
438    SparseCholeskySolver, SparseLuResult, SparseLuSolver, SparseSolver, SymbolicAnalysis,
439};
440
441// Sparse QR factorization (v0.3.0)
442pub mod sparse_qr;
443pub use sparse_qr::{
444    apply_q, apply_qt, extract_q_dense, numerical_rank, sparse_least_squares,
445    sparse_qr as sparse_qr_factorize, SparseLeastSquaresResult, SparseQrConfig, SparseQrResult,
446};
447
448// Sparse matrix reordering algorithms (v0.4.0)
449pub mod reorder;
450pub use reorder::{
451    amd, amd_simple, apply_permutation as reorder_apply_permutation, apply_permutation_csr_array,
452    bandwidth as reorder_bandwidth, cuthill_mckee, cuthill_mckee_full, distance2_color,
453    dsatur_color, greedy_color, nested_dissection, nested_dissection_full,
454    nested_dissection_with_config, profile as reorder_profile,
455    reverse_cuthill_mckee as reorder_rcm, reverse_cuthill_mckee_full, verify_coloring,
456    verify_distance2_coloring, AdjacencyGraph, AmdResult, ColoringResult, CuthillMcKeeResult,
457    GreedyOrdering, NestedDissectionConfig, NestedDissectionResult,
458};
459
460// Unified sparse eigenvalue interface (v0.3.0)
461pub mod sparse_eigen;
462pub use sparse_eigen::{
463    cayley_transform_matvec, check_eigenpairs, compute_residuals, shift_invert_eig, sparse_eig,
464    sparse_eigs, sparse_eigsh, EigenMethod, EigenvalueTarget, SparseEigenConfig, SparseEigenResult,
465    SpectralTransform,
466};
467
468// Quantum-inspired sparse matrix operations (Advanced mode)
469pub mod quantum_inspired_sparse;
470pub use quantum_inspired_sparse::{
471    QuantumProcessorStats, QuantumSparseConfig, QuantumSparseProcessor, QuantumStrategy,
472};
473
474// Neural-adaptive sparse matrix operations (Advanced mode)
475pub mod neural_adaptive_sparse;
476pub use neural_adaptive_sparse::{
477    NeuralAdaptiveConfig, NeuralAdaptiveSparseProcessor, NeuralProcessorStats, OptimizationStrategy,
478};
479
480// Quantum-Neural hybrid optimization (Advanced mode)
481pub mod quantum_neural_hybrid;
482pub use quantum_neural_hybrid::{
483    HybridStrategy, QuantumNeuralConfig, QuantumNeuralHybridProcessor, QuantumNeuralHybridStats,
484};
485
486// Adaptive memory compression for advanced-large sparse matrices (Advanced mode)
487pub mod adaptive_memory_compression;
488pub use adaptive_memory_compression::{
489    AdaptiveCompressionConfig, AdaptiveMemoryCompressor, CompressedMatrix, CompressionAlgorithm,
490    MemoryStats,
491};
492
493// Real-time performance monitoring and adaptation (Advanced mode)
494pub mod realtime_performance_monitor;
495pub use realtime_performance_monitor::{
496    Alert, AlertSeverity, PerformanceMonitorConfig, PerformanceSample, ProcessorType,
497    RealTimePerformanceMonitor,
498};
499
500// Cholesky update/downdate
501pub mod cholesky_update;
502// Distributed sparse operations (SpMV halo exchange, dist AMG)
503pub mod distributed;
504// Learned multigrid smoother
505pub mod learned_smoother;
506// Low-rank sparse update
507pub mod low_rank_update;
508// ML-based preconditioner
509pub mod ml_preconditioner;
510// Parallel AMG
511pub mod parallel_amg;
512
513// Compressed sparse graph algorithms
514pub mod csgraph;
515pub use csgraph::{
516    all_pairs_shortest_path,
517    bellman_ford_single_source,
518    // Centrality measures (v0.2.0)
519    betweenness_centrality,
520    bfs_distances,
521    // Traversal algorithms
522    breadth_first_search,
523    closeness_centrality,
524    // Community detection (v0.2.0)
525    community_detection,
526    compute_laplacianmatrix,
527    connected_components,
528    degree_matrix,
529    depth_first_search,
530    dijkstra_single_source,
531    // Max flow algorithms (v0.2.0)
532    dinic,
533    edmonds_karp,
534    eigenvector_centrality,
535    floyd_warshall,
536    ford_fulkerson,
537    has_path,
538    is_connected,
539    is_laplacian,
540    is_spanning_tree,
541    // Minimum spanning trees
542    kruskal_mst,
543    label_propagation,
544    // Laplacian matrices
545    laplacian,
546    largest_component,
547    louvain_communities,
548    min_cut,
549    minimum_spanning_tree,
550    modularity,
551    num_edges,
552    num_vertices,
553    pagerank,
554    prim_mst,
555    reachable_vertices,
556    reconstruct_path,
557    // Graph algorithms
558    shortest_path,
559    // Shortest path algorithms
560    single_source_shortest_path,
561    spanning_tree_weight,
562    strongly_connected_components,
563    to_adjacency_list,
564    topological_sort,
565    traversegraph,
566    // Connected components
567    undirected_connected_components,
568    // Graph utilities
569    validate_graph,
570    weakly_connected_components,
571    LaplacianType,
572    MSTAlgorithm,
573    // Max flow types (v0.2.0)
574    MaxFlowResult,
575    // Enums and types
576    ShortestPathMethod,
577    TraversalOrder,
578};
579
580// Re-export warnings from scipy for compatibility
581pub struct SparseEfficiencyWarning;
582pub struct SparseWarning;
583
584/// Check if an object is a sparse array
585#[allow(dead_code)]
586pub fn is_sparse_array<T>(obj: &dyn SparseArray<T>) -> bool
587where
588    T: scirs2_core::SparseElement + std::ops::Div<Output = T> + PartialOrd + 'static,
589{
590    sparray::is_sparse(obj)
591}
592
593/// Check if an object is a symmetric sparse array
594#[allow(dead_code)]
595pub fn is_sym_sparse_array<T>(obj: &dyn SymSparseArray<T>) -> bool
596where
597    T: scirs2_core::SparseElement
598        + std::ops::Div<Output = T>
599        + scirs2_core::Float
600        + PartialOrd
601        + 'static,
602{
603    obj.is_symmetric()
604}
605
606/// Check if an object is a sparse matrix (legacy API)
607#[allow(dead_code)]
608pub fn is_sparse_matrix(obj: &dyn std::any::Any) -> bool {
609    obj.is::<CsrMatrix<f64>>()
610        || obj.is::<CscMatrix<f64>>()
611        || obj.is::<CooMatrix<f64>>()
612        || obj.is::<DokMatrix<f64>>()
613        || obj.is::<LilMatrix<f64>>()
614        || obj.is::<DiaMatrix<f64>>()
615        || obj.is::<BsrMatrix<f64>>()
616        || obj.is::<SymCsrMatrix<f64>>()
617        || obj.is::<SymCooMatrix<f64>>()
618        || obj.is::<CsrMatrix<f32>>()
619        || obj.is::<CscMatrix<f32>>()
620        || obj.is::<CooMatrix<f32>>()
621        || obj.is::<DokMatrix<f32>>()
622        || obj.is::<LilMatrix<f32>>()
623        || obj.is::<DiaMatrix<f32>>()
624        || obj.is::<BsrMatrix<f32>>()
625        || obj.is::<SymCsrMatrix<f32>>()
626        || obj.is::<SymCooMatrix<f32>>()
627}
628
629#[cfg(test)]
630mod tests {
631    use super::*;
632    use approx::assert_relative_eq;
633
634    #[test]
635    fn test_csr_array() {
636        let rows = vec![0, 0, 1, 2, 2];
637        let cols = vec![0, 2, 2, 0, 1];
638        let data = vec![1.0, 2.0, 3.0, 4.0, 5.0];
639        let shape = (3, 3);
640
641        let array =
642            CsrArray::from_triplets(&rows, &cols, &data, shape, false).expect("Operation failed");
643
644        assert_eq!(array.shape(), (3, 3));
645        assert_eq!(array.nnz(), 5);
646        assert!(is_sparse_array(&array));
647    }
648
649    #[test]
650    fn test_coo_array() {
651        let rows = vec![0, 0, 1, 2, 2];
652        let cols = vec![0, 2, 2, 0, 1];
653        let data = vec![1.0, 2.0, 3.0, 4.0, 5.0];
654        let shape = (3, 3);
655
656        let array =
657            CooArray::from_triplets(&rows, &cols, &data, shape, false).expect("Operation failed");
658
659        assert_eq!(array.shape(), (3, 3));
660        assert_eq!(array.nnz(), 5);
661        assert!(is_sparse_array(&array));
662    }
663
664    #[test]
665    fn test_dok_array() {
666        let rows = vec![0, 0, 1, 2, 2];
667        let cols = vec![0, 2, 2, 0, 1];
668        let data = vec![1.0, 2.0, 3.0, 4.0, 5.0];
669        let shape = (3, 3);
670
671        let array = DokArray::from_triplets(&rows, &cols, &data, shape).expect("Operation failed");
672
673        assert_eq!(array.shape(), (3, 3));
674        assert_eq!(array.nnz(), 5);
675        assert!(is_sparse_array(&array));
676
677        // Test setting and getting values
678        let mut array = DokArray::<f64>::new((2, 2));
679        array.set(0, 0, 1.0).expect("Operation failed");
680        array.set(1, 1, 2.0).expect("Operation failed");
681
682        assert_eq!(array.get(0, 0), 1.0);
683        assert_eq!(array.get(0, 1), 0.0);
684        assert_eq!(array.get(1, 1), 2.0);
685
686        // Test removing zeros
687        array.set(0, 0, 0.0).expect("Operation failed");
688        assert_eq!(array.nnz(), 1);
689    }
690
691    #[test]
692    fn test_lil_array() {
693        let rows = vec![0, 0, 1, 2, 2];
694        let cols = vec![0, 2, 2, 0, 1];
695        let data = vec![1.0, 2.0, 3.0, 4.0, 5.0];
696        let shape = (3, 3);
697
698        let array = LilArray::from_triplets(&rows, &cols, &data, shape).expect("Operation failed");
699
700        assert_eq!(array.shape(), (3, 3));
701        assert_eq!(array.nnz(), 5);
702        assert!(is_sparse_array(&array));
703
704        // Test setting and getting values
705        let mut array = LilArray::<f64>::new((2, 2));
706        array.set(0, 0, 1.0).expect("Operation failed");
707        array.set(1, 1, 2.0).expect("Operation failed");
708
709        assert_eq!(array.get(0, 0), 1.0);
710        assert_eq!(array.get(0, 1), 0.0);
711        assert_eq!(array.get(1, 1), 2.0);
712
713        // Test sorted indices
714        assert!(array.has_sorted_indices());
715
716        // Test removing zeros
717        array.set(0, 0, 0.0).expect("Operation failed");
718        assert_eq!(array.nnz(), 1);
719    }
720
721    #[test]
722    fn test_dia_array() {
723        use scirs2_core::ndarray::Array1;
724
725        // Create a 3x3 diagonal matrix with main diagonal + upper diagonal
726        let data = vec![
727            Array1::from_vec(vec![1.0, 2.0, 3.0]), // Main diagonal
728            Array1::from_vec(vec![4.0, 5.0, 0.0]), // Upper diagonal
729        ];
730        let offsets = vec![0, 1]; // Main diagonal and k=1
731        let shape = (3, 3);
732
733        let array = DiaArray::new(data, offsets, shape).expect("Operation failed");
734
735        assert_eq!(array.shape(), (3, 3));
736        assert_eq!(array.nnz(), 5); // 3 on main diagonal, 2 on upper diagonal
737        assert!(is_sparse_array(&array));
738
739        // Test values
740        assert_eq!(array.get(0, 0), 1.0);
741        assert_eq!(array.get(1, 1), 2.0);
742        assert_eq!(array.get(2, 2), 3.0);
743        assert_eq!(array.get(0, 1), 4.0);
744        assert_eq!(array.get(1, 2), 5.0);
745        assert_eq!(array.get(0, 2), 0.0);
746
747        // Test from_triplets
748        let rows = vec![0, 0, 1, 1, 2];
749        let cols = vec![0, 1, 1, 2, 2];
750        let data_vec = vec![1.0, 4.0, 2.0, 5.0, 3.0];
751
752        let array2 =
753            DiaArray::from_triplets(&rows, &cols, &data_vec, shape).expect("Operation failed");
754
755        // Should have same values
756        assert_eq!(array2.get(0, 0), 1.0);
757        assert_eq!(array2.get(1, 1), 2.0);
758        assert_eq!(array2.get(2, 2), 3.0);
759        assert_eq!(array2.get(0, 1), 4.0);
760        assert_eq!(array2.get(1, 2), 5.0);
761
762        // Test conversion to other formats
763        let csr = array.to_csr().expect("Operation failed");
764        assert_eq!(csr.nnz(), 5);
765        assert_eq!(csr.get(0, 0), 1.0);
766        assert_eq!(csr.get(0, 1), 4.0);
767    }
768
769    #[test]
770    fn test_format_conversions() {
771        let rows = vec![0, 0, 1, 2, 2];
772        let cols = vec![0, 2, 1, 0, 2];
773        let data = vec![1.0, 2.0, 3.0, 4.0, 5.0];
774        let shape = (3, 3);
775
776        // Create a COO array
777        let coo =
778            CooArray::from_triplets(&rows, &cols, &data, shape, false).expect("Operation failed");
779
780        // Convert to CSR
781        let csr = coo.to_csr().expect("Operation failed");
782
783        // Check values are preserved
784        let coo_dense = coo.to_array();
785        let csr_dense = csr.to_array();
786
787        for i in 0..shape.0 {
788            for j in 0..shape.1 {
789                assert_relative_eq!(coo_dense[[i, j]], csr_dense[[i, j]]);
790            }
791        }
792    }
793
794    #[test]
795    fn test_dot_product() {
796        let rows = vec![0, 0, 1, 2, 2];
797        let cols = vec![0, 2, 1, 0, 2];
798        let data = vec![1.0, 2.0, 3.0, 4.0, 5.0];
799        let shape = (3, 3);
800
801        // Create arrays in different formats
802        let coo =
803            CooArray::from_triplets(&rows, &cols, &data, shape, false).expect("Operation failed");
804        let csr =
805            CsrArray::from_triplets(&rows, &cols, &data, shape, false).expect("Operation failed");
806
807        // Compute dot product (matrix multiplication)
808        let coo_result = coo.dot(&coo).expect("Operation failed");
809        let csr_result = csr.dot(&csr).expect("Operation failed");
810
811        // Check results match
812        let coo_dense = coo_result.to_array();
813        let csr_dense = csr_result.to_array();
814
815        for i in 0..shape.0 {
816            for j in 0..shape.1 {
817                assert_relative_eq!(coo_dense[[i, j]], csr_dense[[i, j]], epsilon = 1e-10);
818            }
819        }
820    }
821
822    #[test]
823    fn test_sym_csr_array() {
824        // Create a symmetric matrix
825        let data = vec![2.0, 1.0, 2.0, 3.0, 0.0, 3.0, 1.0];
826        let indices = vec![0, 0, 1, 2, 0, 1, 2];
827        let indptr = vec![0, 1, 3, 7];
828
829        let sym_matrix =
830            SymCsrMatrix::new(data, indptr, indices, (3, 3)).expect("Operation failed");
831        let sym_array = SymCsrArray::new(sym_matrix);
832
833        assert_eq!(sym_array.shape(), (3, 3));
834        assert!(is_sym_sparse_array(&sym_array));
835
836        // Check values
837        assert_eq!(SparseArray::get(&sym_array, 0, 0), 2.0);
838        assert_eq!(SparseArray::get(&sym_array, 0, 1), 1.0);
839        assert_eq!(SparseArray::get(&sym_array, 1, 0), 1.0); // Symmetric element
840        assert_eq!(SparseArray::get(&sym_array, 1, 2), 3.0);
841        assert_eq!(SparseArray::get(&sym_array, 2, 1), 3.0); // Symmetric element
842
843        // Convert to standard CSR
844        let csr = SymSparseArray::to_csr(&sym_array).expect("Operation failed");
845        assert_eq!(csr.nnz(), 10); // Full matrix with symmetric elements
846    }
847
848    #[test]
849    fn test_sym_coo_array() {
850        // Create a symmetric matrix in COO format
851        let data = vec![2.0, 1.0, 2.0, 3.0, 1.0];
852        let rows = vec![0, 1, 1, 2, 2];
853        let cols = vec![0, 0, 1, 1, 2];
854
855        let sym_matrix = SymCooMatrix::new(data, rows, cols, (3, 3)).expect("Operation failed");
856        let sym_array = SymCooArray::new(sym_matrix);
857
858        assert_eq!(sym_array.shape(), (3, 3));
859        assert!(is_sym_sparse_array(&sym_array));
860
861        // Check values
862        assert_eq!(SparseArray::get(&sym_array, 0, 0), 2.0);
863        assert_eq!(SparseArray::get(&sym_array, 0, 1), 1.0);
864        assert_eq!(SparseArray::get(&sym_array, 1, 0), 1.0); // Symmetric element
865        assert_eq!(SparseArray::get(&sym_array, 1, 2), 3.0);
866        assert_eq!(SparseArray::get(&sym_array, 2, 1), 3.0); // Symmetric element
867
868        // Test from_triplets with enforce symmetry
869        // Input is intentionally asymmetric - will be fixed by enforce_symmetric=true
870        let rows2 = vec![0, 0, 1, 1, 2, 1, 0];
871        let cols2 = vec![0, 1, 1, 2, 2, 0, 2];
872        let data2 = vec![2.0, 1.5, 2.0, 3.5, 1.0, 0.5, 0.0];
873
874        let sym_array2 = SymCooArray::from_triplets(&rows2, &cols2, &data2, (3, 3), true)
875            .expect("Operation failed");
876
877        // Should average the asymmetric values
878        assert_eq!(SparseArray::get(&sym_array2, 0, 1), 1.0); // Average of 1.5 and 0.5
879        assert_eq!(SparseArray::get(&sym_array2, 1, 0), 1.0); // Symmetric element
880        assert_eq!(SparseArray::get(&sym_array2, 0, 2), 0.0); // Zero element
881    }
882
883    #[test]
884    fn test_construct_sym_utils() {
885        // Test creating an identity matrix
886        let eye = construct_sym::eye_sym_array::<f64>(3, "csr").expect("Operation failed");
887
888        assert_eq!(eye.shape(), (3, 3));
889        assert_eq!(SparseArray::get(&*eye, 0, 0), 1.0);
890        assert_eq!(SparseArray::get(&*eye, 1, 1), 1.0);
891        assert_eq!(SparseArray::get(&*eye, 2, 2), 1.0);
892        assert_eq!(SparseArray::get(&*eye, 0, 1), 0.0);
893
894        // Test creating a tridiagonal matrix - with coo format since csr had issues
895        let diag = vec![2.0, 2.0, 2.0];
896        let offdiag = vec![1.0, 1.0];
897
898        let tri =
899            construct_sym::tridiagonal_sym_array(&diag, &offdiag, "coo").expect("Operation failed");
900
901        assert_eq!(tri.shape(), (3, 3));
902        assert_eq!(SparseArray::get(&*tri, 0, 0), 2.0); // Main diagonal
903        assert_eq!(SparseArray::get(&*tri, 1, 1), 2.0);
904        assert_eq!(SparseArray::get(&*tri, 2, 2), 2.0);
905        assert_eq!(SparseArray::get(&*tri, 0, 1), 1.0); // Off-diagonal
906        assert_eq!(SparseArray::get(&*tri, 1, 0), 1.0); // Symmetric element
907        assert_eq!(SparseArray::get(&*tri, 1, 2), 1.0);
908        assert_eq!(SparseArray::get(&*tri, 0, 2), 0.0); // Zero element
909
910        // Test creating a banded matrix
911        let diagonals = vec![
912            vec![2.0, 2.0, 2.0, 2.0, 2.0], // Main diagonal
913            vec![1.0, 1.0, 1.0, 1.0],      // First off-diagonal
914            vec![0.5, 0.5, 0.5],           // Second off-diagonal
915        ];
916
917        let band = construct_sym::banded_sym_array(&diagonals, 5, "csr").expect("Operation failed");
918
919        assert_eq!(band.shape(), (5, 5));
920        assert_eq!(SparseArray::get(&*band, 0, 0), 2.0);
921        assert_eq!(SparseArray::get(&*band, 0, 1), 1.0);
922        assert_eq!(SparseArray::get(&*band, 0, 2), 0.5);
923        assert_eq!(SparseArray::get(&*band, 2, 0), 0.5); // Symmetric element
924    }
925
926    #[test]
927    fn test_sym_conversions() {
928        // Create a symmetric matrix
929        // Lower triangular part only
930        let data = vec![2.0, 1.0, 2.0, 3.0, 1.0];
931        let rows = vec![0, 1, 1, 2, 2];
932        let cols = vec![0, 0, 1, 1, 2];
933
934        let sym_coo = SymCooArray::from_triplets(&rows, &cols, &data, (3, 3), true)
935            .expect("Operation failed");
936
937        // Convert to symmetric CSR
938        let sym_csr = sym_coo.to_sym_csr().expect("Operation failed");
939
940        // Check values are preserved
941        for i in 0..3 {
942            for j in 0..3 {
943                assert_eq!(
944                    SparseArray::get(&sym_coo, i, j),
945                    SparseArray::get(&sym_csr, i, j)
946                );
947            }
948        }
949
950        // Convert to standard formats
951        let csr = SymSparseArray::to_csr(&sym_coo).expect("Operation failed");
952        let coo = SymSparseArray::to_coo(&sym_csr).expect("Operation failed");
953
954        // Check full symmetric matrix in standard formats
955        assert_eq!(csr.nnz(), 7); // Accounts for symmetric pairs
956        assert_eq!(coo.nnz(), 7);
957
958        for i in 0..3 {
959            for j in 0..3 {
960                assert_eq!(SparseArray::get(&csr, i, j), SparseArray::get(&coo, i, j));
961                assert_eq!(
962                    SparseArray::get(&csr, i, j),
963                    SparseArray::get(&sym_csr, i, j)
964                );
965            }
966        }
967    }
968}
scirs2_sparse/lib.rs

scirs2_sparse/
lib.rs
