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.1.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.1.0 (December 29, 2025)
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    lanczos,
212    // Decomposition functions
213    lu_decomposition,
214    matmul,
215    matrix_power,
216    multiply,
217    norm,
218    onenormest,
219    // Eigenvalue functions
220    power_iteration,
221    qr_decomposition,
222    sparse_direct_solve,
223    sparse_lstsq,
224    spsolve,
225    svd_truncated,
226    // SVD functions
227    svds,
228    // TFQMR solver
229    tfqmr,
230    ArpackOptions,
231    // Interfaces
232    AsLinearOperator,
233    // Types from iterative
234    BiCGOptions,
235    BiCGSTABOptions,
236    BiCGSTABResult,
237    // Enhanced operator types
238    BoundaryCondition,
239    CGOptions,
240    CGSOptions,
241    CGSResult,
242    CholeskyResult,
243    ConvolutionMode,
244    ConvolutionOperator,
245    // Operator types
246    DiagonalOperator,
247    EigenResult,
248    EigenvalueMethod,
249    EnhancedDiagonalOperator,
250    EnhancedDifferenceOperator,
251    EnhancedOperatorOptions,
252    EnhancedScaledOperator,
253    EnhancedSumOperator,
254    FiniteDifferenceOperator,
255    GCROTOptions,
256    GCROTResult,
257    GMRESOptions,
258    ICOptions,
259    // Preconditioners
260    ILU0Preconditioner,
261    ILUOptions,
262    IdentityOperator,
263    IterationResult,
264    JacobiPreconditioner,
265    // Decomposition types
266    LUResult,
267    LanczosOptions,
268    LinearOperator,
269    // Eigenvalue types
270    PowerIterationOptions,
271    QRResult,
272    SSORPreconditioner,
273    // SVD types
274    SVDOptions,
275    SVDResult,
276    ScaledIdentityOperator,
277    TFQMROptions,
278    TFQMRResult,
279};
280
281// Format conversions
282pub mod convert;
283
284// Construction utilities
285pub mod construct;
286pub mod construct_sym;
287
288// Combining arrays
289pub mod combine;
290pub use combine::{block_diag, bmat, hstack, kron, kronsum, tril, triu, vstack};
291
292// Index dtype handling utilities
293pub mod index_dtype;
294pub use index_dtype::{can_cast_safely, get_index_dtype, safely_cast_index_arrays};
295
296// Optimized operations for symmetric sparse formats
297pub mod sym_ops;
298pub use sym_ops::{
299    sym_coo_matvec, sym_csr_matvec, sym_csr_quadratic_form, sym_csr_rank1_update, sym_csr_trace,
300};
301
302// GPU-accelerated operations
303pub mod gpu;
304pub mod gpu_kernel_execution;
305pub mod gpu_ops;
306pub mod gpu_spmv_implementation;
307pub use gpu_kernel_execution::{
308    calculate_adaptive_workgroup_size, execute_spmv_kernel, execute_symmetric_spmv_kernel,
309    execute_triangular_solve_kernel, GpuKernelConfig, GpuMemoryManager as GpuKernelMemoryManager,
310    GpuPerformanceProfiler, MemoryStrategy,
311};
312pub use gpu_ops::{
313    gpu_sparse_matvec, gpu_sym_sparse_matvec, AdvancedGpuOps, GpuKernelScheduler, GpuMemoryManager,
314    GpuOptions, GpuProfiler, OptimizedGpuOps,
315};
316pub use gpu_spmv_implementation::GpuSpMV;
317
318// Memory-efficient algorithms and patterns
319pub mod memory_efficient;
320pub use memory_efficient::{
321    streaming_sparse_matvec, CacheAwareOps, MemoryPool, MemoryTracker, OutOfCoreProcessor,
322};
323
324// SIMD-accelerated operations
325pub mod simd_ops;
326pub use simd_ops::{
327    simd_csr_matvec, simd_sparse_elementwise, simd_sparse_linear_combination, simd_sparse_matmul,
328    simd_sparse_norm, simd_sparse_scale, simd_sparse_transpose, ElementwiseOp, SimdOptions,
329};
330
331// Parallel vector operations for iterative solvers
332pub mod parallel_vector_ops;
333pub use parallel_vector_ops::{
334    advanced_sparse_matvec_csr, parallel_axpy, parallel_dot, parallel_linear_combination,
335    parallel_norm2, parallel_sparse_matvec_csr, parallel_vector_add, parallel_vector_copy,
336    parallel_vector_scale, parallel_vector_sub, ParallelVectorOptions,
337};
338
339// Quantum-inspired sparse matrix operations (Advanced mode)
340pub mod quantum_inspired_sparse;
341pub use quantum_inspired_sparse::{
342    QuantumProcessorStats, QuantumSparseConfig, QuantumSparseProcessor, QuantumStrategy,
343};
344
345// Neural-adaptive sparse matrix operations (Advanced mode)
346pub mod neural_adaptive_sparse;
347pub use neural_adaptive_sparse::{
348    NeuralAdaptiveConfig, NeuralAdaptiveSparseProcessor, NeuralProcessorStats, OptimizationStrategy,
349};
350
351// Quantum-Neural hybrid optimization (Advanced mode)
352pub mod quantum_neural_hybrid;
353pub use quantum_neural_hybrid::{
354    HybridStrategy, QuantumNeuralConfig, QuantumNeuralHybridProcessor, QuantumNeuralHybridStats,
355};
356
357// Adaptive memory compression for advanced-large sparse matrices (Advanced mode)
358pub mod adaptive_memory_compression;
359pub use adaptive_memory_compression::{
360    AdaptiveCompressionConfig, AdaptiveMemoryCompressor, CompressedMatrix, CompressionAlgorithm,
361    MemoryStats,
362};
363
364// Real-time performance monitoring and adaptation (Advanced mode)
365pub mod realtime_performance_monitor;
366pub use realtime_performance_monitor::{
367    Alert, AlertSeverity, PerformanceMonitorConfig, PerformanceSample, ProcessorType,
368    RealTimePerformanceMonitor,
369};
370
371// Compressed sparse graph algorithms
372pub mod csgraph;
373pub use csgraph::{
374    all_pairs_shortest_path,
375    bellman_ford_single_source,
376    bfs_distances,
377    // Traversal algorithms
378    breadth_first_search,
379    compute_laplacianmatrix,
380    connected_components,
381    degree_matrix,
382    depth_first_search,
383    dijkstra_single_source,
384    floyd_warshall,
385    has_path,
386    is_connected,
387    is_laplacian,
388    is_spanning_tree,
389    // Minimum spanning trees
390    kruskal_mst,
391    // Laplacian matrices
392    laplacian,
393    largest_component,
394    minimum_spanning_tree,
395    num_edges,
396    num_vertices,
397    prim_mst,
398    reachable_vertices,
399    reconstruct_path,
400    // Graph algorithms
401    shortest_path,
402    // Shortest path algorithms
403    single_source_shortest_path,
404    spanning_tree_weight,
405    strongly_connected_components,
406    to_adjacency_list,
407    topological_sort,
408    traversegraph,
409    // Connected components
410    undirected_connected_components,
411    // Graph utilities
412    validate_graph,
413    weakly_connected_components,
414    LaplacianType,
415    MSTAlgorithm,
416    // Enums and types
417    ShortestPathMethod,
418    TraversalOrder,
419};
420
421// Re-export warnings from scipy for compatibility
422pub struct SparseEfficiencyWarning;
423pub struct SparseWarning;
424
425/// Check if an object is a sparse array
426#[allow(dead_code)]
427pub fn is_sparse_array<T>(obj: &dyn SparseArray<T>) -> bool
428where
429    T: scirs2_core::SparseElement + std::ops::Div<Output = T> + PartialOrd + 'static,
430{
431    sparray::is_sparse(obj)
432}
433
434/// Check if an object is a symmetric sparse array
435#[allow(dead_code)]
436pub fn is_sym_sparse_array<T>(obj: &dyn SymSparseArray<T>) -> bool
437where
438    T: scirs2_core::SparseElement
439        + std::ops::Div<Output = T>
440        + scirs2_core::Float
441        + PartialOrd
442        + 'static,
443{
444    obj.is_symmetric()
445}
446
447/// Check if an object is a sparse matrix (legacy API)
448#[allow(dead_code)]
449pub fn is_sparse_matrix(obj: &dyn std::any::Any) -> bool {
450    obj.is::<CsrMatrix<f64>>()
451        || obj.is::<CscMatrix<f64>>()
452        || obj.is::<CooMatrix<f64>>()
453        || obj.is::<DokMatrix<f64>>()
454        || obj.is::<LilMatrix<f64>>()
455        || obj.is::<DiaMatrix<f64>>()
456        || obj.is::<BsrMatrix<f64>>()
457        || obj.is::<SymCsrMatrix<f64>>()
458        || obj.is::<SymCooMatrix<f64>>()
459        || obj.is::<CsrMatrix<f32>>()
460        || obj.is::<CscMatrix<f32>>()
461        || obj.is::<CooMatrix<f32>>()
462        || obj.is::<DokMatrix<f32>>()
463        || obj.is::<LilMatrix<f32>>()
464        || obj.is::<DiaMatrix<f32>>()
465        || obj.is::<BsrMatrix<f32>>()
466        || obj.is::<SymCsrMatrix<f32>>()
467        || obj.is::<SymCooMatrix<f32>>()
468}
469
470#[cfg(test)]
471mod tests {
472    use super::*;
473    use approx::assert_relative_eq;
474
475    #[test]
476    fn test_csr_array() {
477        let rows = vec![0, 0, 1, 2, 2];
478        let cols = vec![0, 2, 2, 0, 1];
479        let data = vec![1.0, 2.0, 3.0, 4.0, 5.0];
480        let shape = (3, 3);
481
482        let array =
483            CsrArray::from_triplets(&rows, &cols, &data, shape, false).expect("Operation failed");
484
485        assert_eq!(array.shape(), (3, 3));
486        assert_eq!(array.nnz(), 5);
487        assert!(is_sparse_array(&array));
488    }
489
490    #[test]
491    fn test_coo_array() {
492        let rows = vec![0, 0, 1, 2, 2];
493        let cols = vec![0, 2, 2, 0, 1];
494        let data = vec![1.0, 2.0, 3.0, 4.0, 5.0];
495        let shape = (3, 3);
496
497        let array =
498            CooArray::from_triplets(&rows, &cols, &data, shape, false).expect("Operation failed");
499
500        assert_eq!(array.shape(), (3, 3));
501        assert_eq!(array.nnz(), 5);
502        assert!(is_sparse_array(&array));
503    }
504
505    #[test]
506    fn test_dok_array() {
507        let rows = vec![0, 0, 1, 2, 2];
508        let cols = vec![0, 2, 2, 0, 1];
509        let data = vec![1.0, 2.0, 3.0, 4.0, 5.0];
510        let shape = (3, 3);
511
512        let array = DokArray::from_triplets(&rows, &cols, &data, shape).expect("Operation failed");
513
514        assert_eq!(array.shape(), (3, 3));
515        assert_eq!(array.nnz(), 5);
516        assert!(is_sparse_array(&array));
517
518        // Test setting and getting values
519        let mut array = DokArray::<f64>::new((2, 2));
520        array.set(0, 0, 1.0).expect("Operation failed");
521        array.set(1, 1, 2.0).expect("Operation failed");
522
523        assert_eq!(array.get(0, 0), 1.0);
524        assert_eq!(array.get(0, 1), 0.0);
525        assert_eq!(array.get(1, 1), 2.0);
526
527        // Test removing zeros
528        array.set(0, 0, 0.0).expect("Operation failed");
529        assert_eq!(array.nnz(), 1);
530    }
531
532    #[test]
533    fn test_lil_array() {
534        let rows = vec![0, 0, 1, 2, 2];
535        let cols = vec![0, 2, 2, 0, 1];
536        let data = vec![1.0, 2.0, 3.0, 4.0, 5.0];
537        let shape = (3, 3);
538
539        let array = LilArray::from_triplets(&rows, &cols, &data, shape).expect("Operation failed");
540
541        assert_eq!(array.shape(), (3, 3));
542        assert_eq!(array.nnz(), 5);
543        assert!(is_sparse_array(&array));
544
545        // Test setting and getting values
546        let mut array = LilArray::<f64>::new((2, 2));
547        array.set(0, 0, 1.0).expect("Operation failed");
548        array.set(1, 1, 2.0).expect("Operation failed");
549
550        assert_eq!(array.get(0, 0), 1.0);
551        assert_eq!(array.get(0, 1), 0.0);
552        assert_eq!(array.get(1, 1), 2.0);
553
554        // Test sorted indices
555        assert!(array.has_sorted_indices());
556
557        // Test removing zeros
558        array.set(0, 0, 0.0).expect("Operation failed");
559        assert_eq!(array.nnz(), 1);
560    }
561
562    #[test]
563    fn test_dia_array() {
564        use scirs2_core::ndarray::Array1;
565
566        // Create a 3x3 diagonal matrix with main diagonal + upper diagonal
567        let data = vec![
568            Array1::from_vec(vec![1.0, 2.0, 3.0]), // Main diagonal
569            Array1::from_vec(vec![4.0, 5.0, 0.0]), // Upper diagonal
570        ];
571        let offsets = vec![0, 1]; // Main diagonal and k=1
572        let shape = (3, 3);
573
574        let array = DiaArray::new(data, offsets, shape).expect("Operation failed");
575
576        assert_eq!(array.shape(), (3, 3));
577        assert_eq!(array.nnz(), 5); // 3 on main diagonal, 2 on upper diagonal
578        assert!(is_sparse_array(&array));
579
580        // Test values
581        assert_eq!(array.get(0, 0), 1.0);
582        assert_eq!(array.get(1, 1), 2.0);
583        assert_eq!(array.get(2, 2), 3.0);
584        assert_eq!(array.get(0, 1), 4.0);
585        assert_eq!(array.get(1, 2), 5.0);
586        assert_eq!(array.get(0, 2), 0.0);
587
588        // Test from_triplets
589        let rows = vec![0, 0, 1, 1, 2];
590        let cols = vec![0, 1, 1, 2, 2];
591        let data_vec = vec![1.0, 4.0, 2.0, 5.0, 3.0];
592
593        let array2 =
594            DiaArray::from_triplets(&rows, &cols, &data_vec, shape).expect("Operation failed");
595
596        // Should have same values
597        assert_eq!(array2.get(0, 0), 1.0);
598        assert_eq!(array2.get(1, 1), 2.0);
599        assert_eq!(array2.get(2, 2), 3.0);
600        assert_eq!(array2.get(0, 1), 4.0);
601        assert_eq!(array2.get(1, 2), 5.0);
602
603        // Test conversion to other formats
604        let csr = array.to_csr().expect("Operation failed");
605        assert_eq!(csr.nnz(), 5);
606        assert_eq!(csr.get(0, 0), 1.0);
607        assert_eq!(csr.get(0, 1), 4.0);
608    }
609
610    #[test]
611    fn test_format_conversions() {
612        let rows = vec![0, 0, 1, 2, 2];
613        let cols = vec![0, 2, 1, 0, 2];
614        let data = vec![1.0, 2.0, 3.0, 4.0, 5.0];
615        let shape = (3, 3);
616
617        // Create a COO array
618        let coo =
619            CooArray::from_triplets(&rows, &cols, &data, shape, false).expect("Operation failed");
620
621        // Convert to CSR
622        let csr = coo.to_csr().expect("Operation failed");
623
624        // Check values are preserved
625        let coo_dense = coo.to_array();
626        let csr_dense = csr.to_array();
627
628        for i in 0..shape.0 {
629            for j in 0..shape.1 {
630                assert_relative_eq!(coo_dense[[i, j]], csr_dense[[i, j]]);
631            }
632        }
633    }
634
635    #[test]
636    fn test_dot_product() {
637        let rows = vec![0, 0, 1, 2, 2];
638        let cols = vec![0, 2, 1, 0, 2];
639        let data = vec![1.0, 2.0, 3.0, 4.0, 5.0];
640        let shape = (3, 3);
641
642        // Create arrays in different formats
643        let coo =
644            CooArray::from_triplets(&rows, &cols, &data, shape, false).expect("Operation failed");
645        let csr =
646            CsrArray::from_triplets(&rows, &cols, &data, shape, false).expect("Operation failed");
647
648        // Compute dot product (matrix multiplication)
649        let coo_result = coo.dot(&coo).expect("Operation failed");
650        let csr_result = csr.dot(&csr).expect("Operation failed");
651
652        // Check results match
653        let coo_dense = coo_result.to_array();
654        let csr_dense = csr_result.to_array();
655
656        for i in 0..shape.0 {
657            for j in 0..shape.1 {
658                assert_relative_eq!(coo_dense[[i, j]], csr_dense[[i, j]], epsilon = 1e-10);
659            }
660        }
661    }
662
663    #[test]
664    fn test_sym_csr_array() {
665        // Create a symmetric matrix
666        let data = vec![2.0, 1.0, 2.0, 3.0, 0.0, 3.0, 1.0];
667        let indices = vec![0, 0, 1, 2, 0, 1, 2];
668        let indptr = vec![0, 1, 3, 7];
669
670        let sym_matrix =
671            SymCsrMatrix::new(data, indptr, indices, (3, 3)).expect("Operation failed");
672        let sym_array = SymCsrArray::new(sym_matrix);
673
674        assert_eq!(sym_array.shape(), (3, 3));
675        assert!(is_sym_sparse_array(&sym_array));
676
677        // Check values
678        assert_eq!(SparseArray::get(&sym_array, 0, 0), 2.0);
679        assert_eq!(SparseArray::get(&sym_array, 0, 1), 1.0);
680        assert_eq!(SparseArray::get(&sym_array, 1, 0), 1.0); // Symmetric element
681        assert_eq!(SparseArray::get(&sym_array, 1, 2), 3.0);
682        assert_eq!(SparseArray::get(&sym_array, 2, 1), 3.0); // Symmetric element
683
684        // Convert to standard CSR
685        let csr = SymSparseArray::to_csr(&sym_array).expect("Operation failed");
686        assert_eq!(csr.nnz(), 10); // Full matrix with symmetric elements
687    }
688
689    #[test]
690    fn test_sym_coo_array() {
691        // Create a symmetric matrix in COO format
692        let data = vec![2.0, 1.0, 2.0, 3.0, 1.0];
693        let rows = vec![0, 1, 1, 2, 2];
694        let cols = vec![0, 0, 1, 1, 2];
695
696        let sym_matrix = SymCooMatrix::new(data, rows, cols, (3, 3)).expect("Operation failed");
697        let sym_array = SymCooArray::new(sym_matrix);
698
699        assert_eq!(sym_array.shape(), (3, 3));
700        assert!(is_sym_sparse_array(&sym_array));
701
702        // Check values
703        assert_eq!(SparseArray::get(&sym_array, 0, 0), 2.0);
704        assert_eq!(SparseArray::get(&sym_array, 0, 1), 1.0);
705        assert_eq!(SparseArray::get(&sym_array, 1, 0), 1.0); // Symmetric element
706        assert_eq!(SparseArray::get(&sym_array, 1, 2), 3.0);
707        assert_eq!(SparseArray::get(&sym_array, 2, 1), 3.0); // Symmetric element
708
709        // Test from_triplets with enforce symmetry
710        // Input is intentionally asymmetric - will be fixed by enforce_symmetric=true
711        let rows2 = vec![0, 0, 1, 1, 2, 1, 0];
712        let cols2 = vec![0, 1, 1, 2, 2, 0, 2];
713        let data2 = vec![2.0, 1.5, 2.0, 3.5, 1.0, 0.5, 0.0];
714
715        let sym_array2 = SymCooArray::from_triplets(&rows2, &cols2, &data2, (3, 3), true)
716            .expect("Operation failed");
717
718        // Should average the asymmetric values
719        assert_eq!(SparseArray::get(&sym_array2, 0, 1), 1.0); // Average of 1.5 and 0.5
720        assert_eq!(SparseArray::get(&sym_array2, 1, 0), 1.0); // Symmetric element
721        assert_eq!(SparseArray::get(&sym_array2, 0, 2), 0.0); // Zero element
722    }
723
724    #[test]
725    fn test_construct_sym_utils() {
726        // Test creating an identity matrix
727        let eye = construct_sym::eye_sym_array::<f64>(3, "csr").expect("Operation failed");
728
729        assert_eq!(eye.shape(), (3, 3));
730        assert_eq!(SparseArray::get(&*eye, 0, 0), 1.0);
731        assert_eq!(SparseArray::get(&*eye, 1, 1), 1.0);
732        assert_eq!(SparseArray::get(&*eye, 2, 2), 1.0);
733        assert_eq!(SparseArray::get(&*eye, 0, 1), 0.0);
734
735        // Test creating a tridiagonal matrix - with coo format since csr had issues
736        let diag = vec![2.0, 2.0, 2.0];
737        let offdiag = vec![1.0, 1.0];
738
739        let tri =
740            construct_sym::tridiagonal_sym_array(&diag, &offdiag, "coo").expect("Operation failed");
741
742        assert_eq!(tri.shape(), (3, 3));
743        assert_eq!(SparseArray::get(&*tri, 0, 0), 2.0); // Main diagonal
744        assert_eq!(SparseArray::get(&*tri, 1, 1), 2.0);
745        assert_eq!(SparseArray::get(&*tri, 2, 2), 2.0);
746        assert_eq!(SparseArray::get(&*tri, 0, 1), 1.0); // Off-diagonal
747        assert_eq!(SparseArray::get(&*tri, 1, 0), 1.0); // Symmetric element
748        assert_eq!(SparseArray::get(&*tri, 1, 2), 1.0);
749        assert_eq!(SparseArray::get(&*tri, 0, 2), 0.0); // Zero element
750
751        // Test creating a banded matrix
752        let diagonals = vec![
753            vec![2.0, 2.0, 2.0, 2.0, 2.0], // Main diagonal
754            vec![1.0, 1.0, 1.0, 1.0],      // First off-diagonal
755            vec![0.5, 0.5, 0.5],           // Second off-diagonal
756        ];
757
758        let band = construct_sym::banded_sym_array(&diagonals, 5, "csr").expect("Operation failed");
759
760        assert_eq!(band.shape(), (5, 5));
761        assert_eq!(SparseArray::get(&*band, 0, 0), 2.0);
762        assert_eq!(SparseArray::get(&*band, 0, 1), 1.0);
763        assert_eq!(SparseArray::get(&*band, 0, 2), 0.5);
764        assert_eq!(SparseArray::get(&*band, 2, 0), 0.5); // Symmetric element
765    }
766
767    #[test]
768    fn test_sym_conversions() {
769        // Create a symmetric matrix
770        // Lower triangular part only
771        let data = vec![2.0, 1.0, 2.0, 3.0, 1.0];
772        let rows = vec![0, 1, 1, 2, 2];
773        let cols = vec![0, 0, 1, 1, 2];
774
775        let sym_coo = SymCooArray::from_triplets(&rows, &cols, &data, (3, 3), true)
776            .expect("Operation failed");
777
778        // Convert to symmetric CSR
779        let sym_csr = sym_coo.to_sym_csr().expect("Operation failed");
780
781        // Check values are preserved
782        for i in 0..3 {
783            for j in 0..3 {
784                assert_eq!(
785                    SparseArray::get(&sym_coo, i, j),
786                    SparseArray::get(&sym_csr, i, j)
787                );
788            }
789        }
790
791        // Convert to standard formats
792        let csr = SymSparseArray::to_csr(&sym_coo).expect("Operation failed");
793        let coo = SymSparseArray::to_coo(&sym_csr).expect("Operation failed");
794
795        // Check full symmetric matrix in standard formats
796        assert_eq!(csr.nnz(), 7); // Accounts for symmetric pairs
797        assert_eq!(coo.nnz(), 7);
798
799        for i in 0..3 {
800            for j in 0..3 {
801                assert_eq!(SparseArray::get(&csr, i, j), SparseArray::get(&coo, i, j));
802                assert_eq!(
803                    SparseArray::get(&csr, i, j),
804                    SparseArray::get(&sym_csr, i, j)
805                );
806            }
807        }
808    }
809}