Sparse Matrix Library

Sparse Matrix

Sparse Matrix Library

A comprehensive C++17 header-only library for solving systems of linear equations with sparse matrices using advanced iterative methods. The library provides efficient implementations of Krylov subspace methods with support for preconditioning and parallel execution.

Mathematical Foundation

Sparse matrices are matrices in which most elements are zero. This property allows for significant optimization in both storage and computation. The library focuses on iterative methods that are particularly effective for large sparse systems where direct methods would be computationally expensive.

Matrix Representations

The library supports two primary sparse matrix formats:

• Coordinate Format (COO): Stores non-zero elements as triplets (row, column, value)
• Compressed Sparse Row (CSR): Optimized format for efficient row-wise operations

Iterative Methods

The library implements Krylov subspace methods, derivatives of the BiConjugate Gradient method, capable of handling positive definite, negative definite, and indefinite matrices.

Conjugate Gradient Method

The classic CG method, proven to converge for every symmetric positive definite (SPD) matrix. The theoretical guarantee is that the method finds the exact solution in no more than m iterations, where m is the matrix size.

SMM::CSRMatrix matrix;
float* rhs = initializeRHS();
float* initial = initialGuess();
float* result = new float[matrix.getDenseRowCount()];

SMM::SolverStatus status = SMM::ConjugateGradient(
    matrix, rhs, initial, result, 100, 1e-6f
);

BiConjugate Gradient Symmetric

A variant optimized for symmetric matrices. For SPD matrices, it yields identical results to the standard CG method, providing an alternative implementation with different numerical properties.

SMM::SolverStatus status = SMM::BiCGSymmetric(
    matrix, rhs, result, 100, 1e-6f
);

Advanced Methods

• Conjugate Gradient Squared: Transpose-free variant with faster convergence but increased susceptibility to rounding errors
• BiConjugate Gradient Stabilized: Robust transpose-free method with smooth convergence behavior

Preconditioning

Preconditioning dramatically improves convergence rates by transforming the system into an equivalent one with better numerical properties.

Symmetric Gauss-Seidel Preconditioner

A static preconditioner with the form $M = (D + L)^{-1}D(D + U)$, where $D$, $L$, and $U$ represent the diagonal, lower triangular, and upper triangular portions of the matrix respectively.

auto preconditioner = matrix.getPreconditioner(
    SMM::SolverPreconditioner::SYMMETRIC_GAUSS_SEIDEL
);
SMM::SolverStatus status = SMM::BiCGStab(
    matrix, rhs, result, 100, 1e-6f, preconditioner
);

Parallel Implementation

The library supports parallel execution using Intel Threading Building Blocks (TBB), providing significant performance improvements for large-scale problems.

Performance Characteristics

• Scalability: Linear performance improvement with core count for appropriately sized problems
• Memory Efficiency: Minimal overhead for parallel execution
• Load Balancing: Automatic work distribution across available threads

Dependencies and Build System

The library uses modern CMake with FetchContent for dependency management:

• Doctest v2.4.8: Unit testing framework
• TBB v2021.5.0: Parallel execution support (optional)

Build Configuration

mkdir build
cmake -B "./build" -DCMAKE_BUILD_TYPE=Release -DSMM_WITH_TESTS=ON
cd build
cmake --build . --config Release
ctest

CMake Options

• SMM_WITH_MULTITHREADING: Enable TBB-based parallelization (default: ON)
• SMM_WITH_TESTS: Build unit tests (default: OFF)
• SMM_WITH_INSTALL: Generate install target (default: OFF)

Matrix Market Support

Limited support for loading Matrix Market format files, enabling interoperability with scientific computing tools and benchmark datasets.

Applications

This library is particularly useful in computational domains where large sparse systems arise naturally.

Use Cases

• Finite Element Analysis: Structural engineering, fluid dynamics
• Graph Algorithms: Network analysis, PageRank computations
• Machine Learning: Large-scale optimization problems
• Scientific Computing: Differential equation discretizations

Implementation Notes

The header-only design simplifies integration while maintaining high performance through template specialization and compiler optimizations.