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.2"
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.4.2 (March 27, 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_preconditioner;
352pub mod gpu_spmv_implementation;
353pub use gpu_kernel_execution::{
354    calculate_adaptive_workgroup_size, execute_spmv_kernel, execute_symmetric_spmv_kernel,
355    execute_triangular_solve_kernel, GpuKernelConfig, GpuMemoryManager as GpuKernelMemoryManager,
356    GpuPerformanceProfiler, MemoryStrategy,
357};
358pub use gpu_ops::{
359    gpu_sparse_matvec, gpu_sym_sparse_matvec, AdvancedGpuOps, GpuKernelScheduler, GpuMemoryManager,
360    GpuOptions, GpuProfiler, OptimizedGpuOps,
361};
362pub use gpu_spmv_implementation::GpuSpMV;
363
364// Memory-efficient algorithms and patterns
365pub mod memory_efficient;
366pub use memory_efficient::{
367    streaming_sparse_matvec, CacheAwareOps, MemoryPool, MemoryTracker, OutOfCoreProcessor,
368};
369
370// SIMD-accelerated operations
371pub mod simd_ops;
372pub use simd_ops::{
373    simd_csr_matvec, simd_sparse_elementwise, simd_sparse_linear_combination, simd_sparse_matmul,
374    simd_sparse_norm, simd_sparse_scale, simd_sparse_transpose, ElementwiseOp, SimdOptions,
375};
376
377// Parallel vector operations for iterative solvers
378pub mod parallel_vector_ops;
379pub use parallel_vector_ops::{
380    advanced_sparse_matvec_csr, parallel_axpy, parallel_dot, parallel_linear_combination,
381    parallel_norm2, parallel_sparse_matvec_csr, parallel_vector_add, parallel_vector_copy,
382    parallel_vector_scale, parallel_vector_sub, ParallelVectorOptions,
383};
384
385// Enhanced iterative solvers with preconditioners and sparse utilities (v0.3.0)
386pub mod iterative_solvers;
387pub use iterative_solvers::{
388    // Solvers
389    bicgstab as enhanced_bicgstab,
390    cg as enhanced_cg,
391    chebyshev,
392    // Sparse utility functions
393    estimate_spectral_radius,
394    gmres as enhanced_gmres,
395    sparse_diagonal,
396    sparse_norm,
397    sparse_trace,
398    // Preconditioners (Array1-based interface)
399    ILU0Preconditioner as EnhancedILU0Preconditioner,
400    // Configuration and result types
401    IterativeSolverConfig,
402    JacobiPreconditioner as EnhancedJacobiPreconditioner,
403    NormType,
404    Preconditioner,
405    SSORPreconditioner as EnhancedSSORPreconditioner,
406    SolverResult,
407};
408
409// LOBPCG eigensolver (v0.3.0)
410pub mod lobpcg;
411pub use lobpcg::{
412    lobpcg as lobpcg_eigensolver, lobpcg_generalised, EigenTarget, LobpcgConfig, LobpcgResult,
413};
414
415// Advanced Krylov subspace eigensolvers (v0.3.0)
416pub mod krylov;
417pub use krylov::{
418    iram as krylov_iram, thick_restart_lanczos, IramConfig, KrylovEigenResult,
419    ThickRestartLanczosConfig, WhichEigenvalues,
420};
421
422// Sparse matrix utilities (v0.3.0)
423pub mod sparse_utils;
424pub use sparse_utils::{
425    condest_1norm, permute_matrix, reverse_cuthill_mckee, sparse_add, sparse_extract_diagonal,
426    sparse_identity, sparse_kronecker, sparse_matrix_norm, sparse_matrix_trace, sparse_scale,
427    sparse_sub, sparse_transpose, spgemm, RcmResult, SparseNorm,
428};
429
430// Incomplete factorization preconditioners (v0.3.0)
431pub mod incomplete_factorizations;
432pub use incomplete_factorizations::{Ic0, Ilu0 as Ilu0Enhanced, IluK, Ilut, IlutConfig, Milu};
433
434// Sparse direct solvers (v0.3.0)
435pub mod direct_solver;
436pub use direct_solver::{
437    amd_ordering, elimination_tree, inverse_perm, nested_dissection_ordering,
438    sparse_cholesky_solve, sparse_lu_solve, symbolic_cholesky, SparseCholResult,
439    SparseCholeskySolver, SparseLuResult, SparseLuSolver, SparseSolver, SymbolicAnalysis,
440};
441
442// Sparse QR factorization (v0.3.0)
443pub mod sparse_qr;
444pub use sparse_qr::{
445    apply_q, apply_qt, extract_q_dense, numerical_rank, sparse_least_squares,
446    sparse_qr as sparse_qr_factorize, SparseLeastSquaresResult, SparseQrConfig, SparseQrResult,
447};
448
449// Sparse matrix reordering algorithms (v0.4.0)
450pub mod reorder;
451pub use reorder::{
452    amd, amd_simple, apply_permutation as reorder_apply_permutation, apply_permutation_csr_array,
453    bandwidth as reorder_bandwidth, cuthill_mckee, cuthill_mckee_full, distance2_color,
454    dsatur_color, greedy_color, nested_dissection, nested_dissection_full,
455    nested_dissection_with_config, profile as reorder_profile,
456    reverse_cuthill_mckee as reorder_rcm, reverse_cuthill_mckee_full, verify_coloring,
457    verify_distance2_coloring, AdjacencyGraph, AmdResult, ColoringResult, CuthillMcKeeResult,
458    GreedyOrdering, NestedDissectionConfig, NestedDissectionResult,
459};
460
461// Unified sparse eigenvalue interface (v0.3.0)
462pub mod sparse_eigen;
463pub use sparse_eigen::{
464    cayley_transform_matvec, check_eigenpairs, compute_residuals, shift_invert_eig, sparse_eig,
465    sparse_eigs, sparse_eigsh, EigenMethod, EigenvalueTarget, SparseEigenConfig, SparseEigenResult,
466    SpectralTransform,
467};
468
469// Quantum-inspired sparse matrix operations (Advanced mode)
470pub mod quantum_inspired_sparse;
471pub use quantum_inspired_sparse::{
472    QuantumProcessorStats, QuantumSparseConfig, QuantumSparseProcessor, QuantumStrategy,
473};
474
475// Neural-adaptive sparse matrix operations (Advanced mode)
476pub mod neural_adaptive_sparse;
477pub use neural_adaptive_sparse::{
478    NeuralAdaptiveConfig, NeuralAdaptiveSparseProcessor, NeuralProcessorStats, OptimizationStrategy,
479};
480
481// Quantum-Neural hybrid optimization (Advanced mode)
482pub mod quantum_neural_hybrid;
483pub use quantum_neural_hybrid::{
484    HybridStrategy, QuantumNeuralConfig, QuantumNeuralHybridProcessor, QuantumNeuralHybridStats,
485};
486
487// Adaptive memory compression for advanced-large sparse matrices (Advanced mode)
488pub mod adaptive_memory_compression;
489pub use adaptive_memory_compression::{
490    AdaptiveCompressionConfig, AdaptiveMemoryCompressor, CompressedMatrix, CompressionAlgorithm,
491    MemoryStats,
492};
493
494// Real-time performance monitoring and adaptation (Advanced mode)
495pub mod realtime_performance_monitor;
496pub use realtime_performance_monitor::{
497    Alert, AlertSeverity, PerformanceMonitorConfig, PerformanceSample, ProcessorType,
498    RealTimePerformanceMonitor,
499};
500
501// Cholesky update/downdate
502pub mod cholesky_update;
503// Distributed sparse operations (SpMV halo exchange, dist AMG)
504pub mod distributed;
505// Learned multigrid smoother
506pub mod learned_smoother;
507// Low-rank sparse update
508pub mod low_rank_update;
509// ML-based preconditioner
510pub mod ml_preconditioner;
511// Parallel AMG
512pub mod parallel_amg;
513
514// Compressed sparse graph algorithms
515pub mod csgraph;
516pub use csgraph::{
517    all_pairs_shortest_path,
518    bellman_ford_single_source,
519    // Centrality measures (v0.2.0)
520    betweenness_centrality,
521    bfs_distances,
522    // Traversal algorithms
523    breadth_first_search,
524    closeness_centrality,
525    // Community detection (v0.2.0)
526    community_detection,
527    compute_laplacianmatrix,
528    connected_components,
529    degree_matrix,
530    depth_first_search,
531    dijkstra_single_source,
532    // Max flow algorithms (v0.2.0)
533    dinic,
534    edmonds_karp,
535    eigenvector_centrality,
536    floyd_warshall,
537    ford_fulkerson,
538    has_path,
539    is_connected,
540    is_laplacian,
541    is_spanning_tree,
542    // Minimum spanning trees
543    kruskal_mst,
544    label_propagation,
545    // Laplacian matrices
546    laplacian,
547    largest_component,
548    louvain_communities,
549    min_cut,
550    minimum_spanning_tree,
551    modularity,
552    num_edges,
553    num_vertices,
554    pagerank,
555    prim_mst,
556    reachable_vertices,
557    reconstruct_path,
558    // Graph algorithms
559    shortest_path,
560    // Shortest path algorithms
561    single_source_shortest_path,
562    spanning_tree_weight,
563    strongly_connected_components,
564    to_adjacency_list,
565    topological_sort,
566    traversegraph,
567    // Connected components
568    undirected_connected_components,
569    // Graph utilities
570    validate_graph,
571    weakly_connected_components,
572    LaplacianType,
573    MSTAlgorithm,
574    // Max flow types (v0.2.0)
575    MaxFlowResult,
576    // Enums and types
577    ShortestPathMethod,
578    TraversalOrder,
579};
580
581// Re-export warnings from scipy for compatibility
582pub struct SparseEfficiencyWarning;
583pub struct SparseWarning;
584
585/// Check if an object is a sparse array
586#[allow(dead_code)]
587pub fn is_sparse_array<T>(obj: &dyn SparseArray<T>) -> bool
588where
589    T: scirs2_core::SparseElement + std::ops::Div<Output = T> + PartialOrd + 'static,
590{
591    sparray::is_sparse(obj)
592}
593
594/// Check if an object is a symmetric sparse array
595#[allow(dead_code)]
596pub fn is_sym_sparse_array<T>(obj: &dyn SymSparseArray<T>) -> bool
597where
598    T: scirs2_core::SparseElement
599        + std::ops::Div<Output = T>
600        + scirs2_core::Float
601        + PartialOrd
602        + 'static,
603{
604    obj.is_symmetric()
605}
606
607/// Check if an object is a sparse matrix (legacy API)
608#[allow(dead_code)]
609pub fn is_sparse_matrix(obj: &dyn std::any::Any) -> bool {
610    obj.is::<CsrMatrix<f64>>()
611        || obj.is::<CscMatrix<f64>>()
612        || obj.is::<CooMatrix<f64>>()
613        || obj.is::<DokMatrix<f64>>()
614        || obj.is::<LilMatrix<f64>>()
615        || obj.is::<DiaMatrix<f64>>()
616        || obj.is::<BsrMatrix<f64>>()
617        || obj.is::<SymCsrMatrix<f64>>()
618        || obj.is::<SymCooMatrix<f64>>()
619        || obj.is::<CsrMatrix<f32>>()
620        || obj.is::<CscMatrix<f32>>()
621        || obj.is::<CooMatrix<f32>>()
622        || obj.is::<DokMatrix<f32>>()
623        || obj.is::<LilMatrix<f32>>()
624        || obj.is::<DiaMatrix<f32>>()
625        || obj.is::<BsrMatrix<f32>>()
626        || obj.is::<SymCsrMatrix<f32>>()
627        || obj.is::<SymCooMatrix<f32>>()
628}
629
630#[cfg(test)]
631mod tests {
632    use super::*;
633    use approx::assert_relative_eq;
634
635    #[test]
636    fn test_csr_array() {
637        let rows = vec![0, 0, 1, 2, 2];
638        let cols = vec![0, 2, 2, 0, 1];
639        let data = vec![1.0, 2.0, 3.0, 4.0, 5.0];
640        let shape = (3, 3);
641
642        let array =
643            CsrArray::from_triplets(&rows, &cols, &data, shape, false).expect("Operation failed");
644
645        assert_eq!(array.shape(), (3, 3));
646        assert_eq!(array.nnz(), 5);
647        assert!(is_sparse_array(&array));
648    }
649
650    #[test]
651    fn test_coo_array() {
652        let rows = vec![0, 0, 1, 2, 2];
653        let cols = vec![0, 2, 2, 0, 1];
654        let data = vec![1.0, 2.0, 3.0, 4.0, 5.0];
655        let shape = (3, 3);
656
657        let array =
658            CooArray::from_triplets(&rows, &cols, &data, shape, false).expect("Operation failed");
659
660        assert_eq!(array.shape(), (3, 3));
661        assert_eq!(array.nnz(), 5);
662        assert!(is_sparse_array(&array));
663    }
664
665    #[test]
666    fn test_dok_array() {
667        let rows = vec![0, 0, 1, 2, 2];
668        let cols = vec![0, 2, 2, 0, 1];
669        let data = vec![1.0, 2.0, 3.0, 4.0, 5.0];
670        let shape = (3, 3);
671
672        let array = DokArray::from_triplets(&rows, &cols, &data, shape).expect("Operation failed");
673
674        assert_eq!(array.shape(), (3, 3));
675        assert_eq!(array.nnz(), 5);
676        assert!(is_sparse_array(&array));
677
678        // Test setting and getting values
679        let mut array = DokArray::<f64>::new((2, 2));
680        array.set(0, 0, 1.0).expect("Operation failed");
681        array.set(1, 1, 2.0).expect("Operation failed");
682
683        assert_eq!(array.get(0, 0), 1.0);
684        assert_eq!(array.get(0, 1), 0.0);
685        assert_eq!(array.get(1, 1), 2.0);
686
687        // Test removing zeros
688        array.set(0, 0, 0.0).expect("Operation failed");
689        assert_eq!(array.nnz(), 1);
690    }
691
692    #[test]
693    fn test_lil_array() {
694        let rows = vec![0, 0, 1, 2, 2];
695        let cols = vec![0, 2, 2, 0, 1];
696        let data = vec![1.0, 2.0, 3.0, 4.0, 5.0];
697        let shape = (3, 3);
698
699        let array = LilArray::from_triplets(&rows, &cols, &data, shape).expect("Operation failed");
700
701        assert_eq!(array.shape(), (3, 3));
702        assert_eq!(array.nnz(), 5);
703        assert!(is_sparse_array(&array));
704
705        // Test setting and getting values
706        let mut array = LilArray::<f64>::new((2, 2));
707        array.set(0, 0, 1.0).expect("Operation failed");
708        array.set(1, 1, 2.0).expect("Operation failed");
709
710        assert_eq!(array.get(0, 0), 1.0);
711        assert_eq!(array.get(0, 1), 0.0);
712        assert_eq!(array.get(1, 1), 2.0);
713
714        // Test sorted indices
715        assert!(array.has_sorted_indices());
716
717        // Test removing zeros
718        array.set(0, 0, 0.0).expect("Operation failed");
719        assert_eq!(array.nnz(), 1);
720    }
721
722    #[test]
723    fn test_dia_array() {
724        use scirs2_core::ndarray::Array1;
725
726        // Create a 3x3 diagonal matrix with main diagonal + upper diagonal
727        let data = vec![
728            Array1::from_vec(vec![1.0, 2.0, 3.0]), // Main diagonal
729            Array1::from_vec(vec![4.0, 5.0, 0.0]), // Upper diagonal
730        ];
731        let offsets = vec![0, 1]; // Main diagonal and k=1
732        let shape = (3, 3);
733
734        let array = DiaArray::new(data, offsets, shape).expect("Operation failed");
735
736        assert_eq!(array.shape(), (3, 3));
737        assert_eq!(array.nnz(), 5); // 3 on main diagonal, 2 on upper diagonal
738        assert!(is_sparse_array(&array));
739
740        // Test values
741        assert_eq!(array.get(0, 0), 1.0);
742        assert_eq!(array.get(1, 1), 2.0);
743        assert_eq!(array.get(2, 2), 3.0);
744        assert_eq!(array.get(0, 1), 4.0);
745        assert_eq!(array.get(1, 2), 5.0);
746        assert_eq!(array.get(0, 2), 0.0);
747
748        // Test from_triplets
749        let rows = vec![0, 0, 1, 1, 2];
750        let cols = vec![0, 1, 1, 2, 2];
751        let data_vec = vec![1.0, 4.0, 2.0, 5.0, 3.0];
752
753        let array2 =
754            DiaArray::from_triplets(&rows, &cols, &data_vec, shape).expect("Operation failed");
755
756        // Should have same values
757        assert_eq!(array2.get(0, 0), 1.0);
758        assert_eq!(array2.get(1, 1), 2.0);
759        assert_eq!(array2.get(2, 2), 3.0);
760        assert_eq!(array2.get(0, 1), 4.0);
761        assert_eq!(array2.get(1, 2), 5.0);
762
763        // Test conversion to other formats
764        let csr = array.to_csr().expect("Operation failed");
765        assert_eq!(csr.nnz(), 5);
766        assert_eq!(csr.get(0, 0), 1.0);
767        assert_eq!(csr.get(0, 1), 4.0);
768    }
769
770    #[test]
771    fn test_format_conversions() {
772        let rows = vec![0, 0, 1, 2, 2];
773        let cols = vec![0, 2, 1, 0, 2];
774        let data = vec![1.0, 2.0, 3.0, 4.0, 5.0];
775        let shape = (3, 3);
776
777        // Create a COO array
778        let coo =
779            CooArray::from_triplets(&rows, &cols, &data, shape, false).expect("Operation failed");
780
781        // Convert to CSR
782        let csr = coo.to_csr().expect("Operation failed");
783
784        // Check values are preserved
785        let coo_dense = coo.to_array();
786        let csr_dense = csr.to_array();
787
788        for i in 0..shape.0 {
789            for j in 0..shape.1 {
790                assert_relative_eq!(coo_dense[[i, j]], csr_dense[[i, j]]);
791            }
792        }
793    }
794
795    #[test]
796    fn test_dot_product() {
797        let rows = vec![0, 0, 1, 2, 2];
798        let cols = vec![0, 2, 1, 0, 2];
799        let data = vec![1.0, 2.0, 3.0, 4.0, 5.0];
800        let shape = (3, 3);
801
802        // Create arrays in different formats
803        let coo =
804            CooArray::from_triplets(&rows, &cols, &data, shape, false).expect("Operation failed");
805        let csr =
806            CsrArray::from_triplets(&rows, &cols, &data, shape, false).expect("Operation failed");
807
808        // Compute dot product (matrix multiplication)
809        let coo_result = coo.dot(&coo).expect("Operation failed");
810        let csr_result = csr.dot(&csr).expect("Operation failed");
811
812        // Check results match
813        let coo_dense = coo_result.to_array();
814        let csr_dense = csr_result.to_array();
815
816        for i in 0..shape.0 {
817            for j in 0..shape.1 {
818                assert_relative_eq!(coo_dense[[i, j]], csr_dense[[i, j]], epsilon = 1e-10);
819            }
820        }
821    }
822
823    #[test]
824    fn test_sym_csr_array() {
825        // Create a symmetric matrix
826        let data = vec![2.0, 1.0, 2.0, 3.0, 0.0, 3.0, 1.0];
827        let indices = vec![0, 0, 1, 2, 0, 1, 2];
828        let indptr = vec![0, 1, 3, 7];
829
830        let sym_matrix =
831            SymCsrMatrix::new(data, indptr, indices, (3, 3)).expect("Operation failed");
832        let sym_array = SymCsrArray::new(sym_matrix);
833
834        assert_eq!(sym_array.shape(), (3, 3));
835        assert!(is_sym_sparse_array(&sym_array));
836
837        // Check values
838        assert_eq!(SparseArray::get(&sym_array, 0, 0), 2.0);
839        assert_eq!(SparseArray::get(&sym_array, 0, 1), 1.0);
840        assert_eq!(SparseArray::get(&sym_array, 1, 0), 1.0); // Symmetric element
841        assert_eq!(SparseArray::get(&sym_array, 1, 2), 3.0);
842        assert_eq!(SparseArray::get(&sym_array, 2, 1), 3.0); // Symmetric element
843
844        // Convert to standard CSR
845        let csr = SymSparseArray::to_csr(&sym_array).expect("Operation failed");
846        assert_eq!(csr.nnz(), 10); // Full matrix with symmetric elements
847    }
848
849    #[test]
850    fn test_sym_coo_array() {
851        // Create a symmetric matrix in COO format
852        let data = vec![2.0, 1.0, 2.0, 3.0, 1.0];
853        let rows = vec![0, 1, 1, 2, 2];
854        let cols = vec![0, 0, 1, 1, 2];
855
856        let sym_matrix = SymCooMatrix::new(data, rows, cols, (3, 3)).expect("Operation failed");
857        let sym_array = SymCooArray::new(sym_matrix);
858
859        assert_eq!(sym_array.shape(), (3, 3));
860        assert!(is_sym_sparse_array(&sym_array));
861
862        // Check values
863        assert_eq!(SparseArray::get(&sym_array, 0, 0), 2.0);
864        assert_eq!(SparseArray::get(&sym_array, 0, 1), 1.0);
865        assert_eq!(SparseArray::get(&sym_array, 1, 0), 1.0); // Symmetric element
866        assert_eq!(SparseArray::get(&sym_array, 1, 2), 3.0);
867        assert_eq!(SparseArray::get(&sym_array, 2, 1), 3.0); // Symmetric element
868
869        // Test from_triplets with enforce symmetry
870        // Input is intentionally asymmetric - will be fixed by enforce_symmetric=true
871        let rows2 = vec![0, 0, 1, 1, 2, 1, 0];
872        let cols2 = vec![0, 1, 1, 2, 2, 0, 2];
873        let data2 = vec![2.0, 1.5, 2.0, 3.5, 1.0, 0.5, 0.0];
874
875        let sym_array2 = SymCooArray::from_triplets(&rows2, &cols2, &data2, (3, 3), true)
876            .expect("Operation failed");
877
878        // Should average the asymmetric values
879        assert_eq!(SparseArray::get(&sym_array2, 0, 1), 1.0); // Average of 1.5 and 0.5
880        assert_eq!(SparseArray::get(&sym_array2, 1, 0), 1.0); // Symmetric element
881        assert_eq!(SparseArray::get(&sym_array2, 0, 2), 0.0); // Zero element
882    }
883
884    #[test]
885    fn test_construct_sym_utils() {
886        // Test creating an identity matrix
887        let eye = construct_sym::eye_sym_array::<f64>(3, "csr").expect("Operation failed");
888
889        assert_eq!(eye.shape(), (3, 3));
890        assert_eq!(SparseArray::get(&*eye, 0, 0), 1.0);
891        assert_eq!(SparseArray::get(&*eye, 1, 1), 1.0);
892        assert_eq!(SparseArray::get(&*eye, 2, 2), 1.0);
893        assert_eq!(SparseArray::get(&*eye, 0, 1), 0.0);
894
895        // Test creating a tridiagonal matrix - with coo format since csr had issues
896        let diag = vec![2.0, 2.0, 2.0];
897        let offdiag = vec![1.0, 1.0];
898
899        let tri =
900            construct_sym::tridiagonal_sym_array(&diag, &offdiag, "coo").expect("Operation failed");
901
902        assert_eq!(tri.shape(), (3, 3));
903        assert_eq!(SparseArray::get(&*tri, 0, 0), 2.0); // Main diagonal
904        assert_eq!(SparseArray::get(&*tri, 1, 1), 2.0);
905        assert_eq!(SparseArray::get(&*tri, 2, 2), 2.0);
906        assert_eq!(SparseArray::get(&*tri, 0, 1), 1.0); // Off-diagonal
907        assert_eq!(SparseArray::get(&*tri, 1, 0), 1.0); // Symmetric element
908        assert_eq!(SparseArray::get(&*tri, 1, 2), 1.0);
909        assert_eq!(SparseArray::get(&*tri, 0, 2), 0.0); // Zero element
910
911        // Test creating a banded matrix
912        let diagonals = vec![
913            vec![2.0, 2.0, 2.0, 2.0, 2.0], // Main diagonal
914            vec![1.0, 1.0, 1.0, 1.0],      // First off-diagonal
915            vec![0.5, 0.5, 0.5],           // Second off-diagonal
916        ];
917
918        let band = construct_sym::banded_sym_array(&diagonals, 5, "csr").expect("Operation failed");
919
920        assert_eq!(band.shape(), (5, 5));
921        assert_eq!(SparseArray::get(&*band, 0, 0), 2.0);
922        assert_eq!(SparseArray::get(&*band, 0, 1), 1.0);
923        assert_eq!(SparseArray::get(&*band, 0, 2), 0.5);
924        assert_eq!(SparseArray::get(&*band, 2, 0), 0.5); // Symmetric element
925    }
926
927    #[test]
928    fn test_sym_conversions() {
929        // Create a symmetric matrix
930        // Lower triangular part only
931        let data = vec![2.0, 1.0, 2.0, 3.0, 1.0];
932        let rows = vec![0, 1, 1, 2, 2];
933        let cols = vec![0, 0, 1, 1, 2];
934
935        let sym_coo = SymCooArray::from_triplets(&rows, &cols, &data, (3, 3), true)
936            .expect("Operation failed");
937
938        // Convert to symmetric CSR
939        let sym_csr = sym_coo.to_sym_csr().expect("Operation failed");
940
941        // Check values are preserved
942        for i in 0..3 {
943            for j in 0..3 {
944                assert_eq!(
945                    SparseArray::get(&sym_coo, i, j),
946                    SparseArray::get(&sym_csr, i, j)
947                );
948            }
949        }
950
951        // Convert to standard formats
952        let csr = SymSparseArray::to_csr(&sym_coo).expect("Operation failed");
953        let coo = SymSparseArray::to_coo(&sym_csr).expect("Operation failed");
954
955        // Check full symmetric matrix in standard formats
956        assert_eq!(csr.nnz(), 7); // Accounts for symmetric pairs
957        assert_eq!(coo.nnz(), 7);
958
959        for i in 0..3 {
960            for j in 0..3 {
961                assert_eq!(SparseArray::get(&csr, i, j), SparseArray::get(&coo, i, j));
962                assert_eq!(
963                    SparseArray::get(&csr, i, j),
964                    SparseArray::get(&sym_csr, i, j)
965                );
966            }
967        }
968    }
969}
scirs2_sparse/lib.rs

scirs2_sparse/
lib.rs
