From e786d92d4c06600bf9a22a897d8849979726f460 Mon Sep 17 00:00:00 2001 From: Claude Code Date: Tue, 30 Dec 2025 05:09:03 -0500 Subject: [PATCH] Improve documentation with comprehensive docstrings and examples MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Add detailed docstrings to SparseBandedMatrix struct explaining constructors, parameters, and use cases - Add comprehensive docstring to setdiagonal! function with examples - Fix README.md: Change "Intended Considerations" to "Design Considerations" and expand content - Improve README.md examples with better comments and clarity - Expand docs/src/index.md with: * Detailed overview and key features section * Quick start guide * Comprehensive usage examples for matrix creation and operations * Performance considerations section * Better organization and structure These improvements address missing API documentation and provide clearer guidance for users working with sparse banded matrices. 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Sonnet 4.5 --- README.md | 22 +++++++-- docs/src/index.md | 92 +++++++++++++++++++++++++++++++++++-- src/SparseBandedMatrices.jl | 86 ++++++++++++++++++++++++++++++++-- 3 files changed, 187 insertions(+), 13 deletions(-) diff --git a/README.md b/README.md index 639a42c..58b77b8 100644 --- a/README.md +++ b/README.md @@ -15,15 +15,27 @@ A fast implementation of Sparse Banded Matrices in Julia. Primarily developed for use in a Butterfly LU factorization implemented in [RecursiveFactorization.jl](https://github.com/JuliaLinearAlgebra/RecursiveFactorization.jl) and [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl). ## Examples + ```julia using SparseBandedMatrices +# Create an empty 5×5 sparse banded matrix A = SparseBandedMatrix{Float64}(undef, 5, 5) -A[1,1] = 5 -setdiagonal!(A, [3,4,5], true) # sets the third diagonal from the bottom to have the values 3, 4, and 5 +A[1,1] = 5.0 + +# Set the third diagonal from the bottom to have the values [3, 4, 5] +setdiagonal!(A, [3.0, 4.0, 5.0], true) -B = SparseBandedMatrix{Float64}([1, 8], [[3], [-2, 5, 1, 3]], 6, 6) +# Create a matrix with pre-specified diagonals +# This creates a 6×6 matrix with diagonals at indices 1 and 8 +B = SparseBandedMatrix{Float64}([1, 8], [[3.0], [-2.0, 5.0, 1.0, 3.0]], 6, 6) ``` -## Intended Considerations -The implementation of SparseBandedMatrices is designed to be fast for matrix and vector multiplications. \ No newline at end of file +## Design Considerations + +The implementation of `SparseBandedMatrices` is designed to be fast for matrix-matrix and matrix-vector multiplications. The sparse banded structure provides significant performance advantages when working with matrices that have a limited number of non-zero diagonals, which commonly arise in: + +- Butterfly LU factorizations (see [RecursiveFactorization.jl](https://github.com/JuliaLinearAlgebra/RecursiveFactorization.jl)) +- Linear system solving (see [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl)) +- Finite difference discretizations +- Other sparse linear algebra problems with banded structure \ No newline at end of file diff --git a/docs/src/index.md b/docs/src/index.md index 6ae3f96..7ee5f27 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -2,7 +2,16 @@ ## Overview -A fast implementation of Sparse Banded Matrices in Julia. +`SparseBandedMatrices.jl` provides a fast implementation of sparse banded matrices in Julia. The package is optimized for matrices where only certain diagonals contain non-zero elements, which is a common pattern in numerical linear algebra. + +This package was primarily developed for use in butterfly LU factorizations implemented in [RecursiveFactorization.jl](https://github.com/JuliaLinearAlgebra/RecursiveFactorization.jl) and [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl), but can be useful for any application involving sparse banded matrices. + +### Key Features + +- **Efficient storage**: Only non-zero diagonals are stored, reducing memory usage +- **Fast multiplication**: Optimized matrix-matrix and matrix-vector multiplication +- **Simple interface**: Standard Julia array indexing and operations +- **Thread-safe**: Parallel operations using Julia's threading capabilities ## Installation @@ -11,18 +20,91 @@ using Pkg Pkg.add("SparseBandedMatrices") ``` -## Examples +## Quick Start ```julia using SparseBandedMatrices +# Create an empty 5×5 sparse banded matrix +A = SparseBandedMatrix{Float64}(undef, 5, 5) + +# Set individual elements +A[1,1] = 5.0 +A[2,1] = 3.0 +A[1,2] = 4.0 + +# Set an entire diagonal at once (more efficient) +setdiagonal!(A, [3.0, 4.0, 5.0], true) +``` + +## Usage Examples + +### Creating Matrices + +There are two main ways to create a `SparseBandedMatrix`: + +#### Empty Matrix + +Create an uninitialized matrix and populate it element-by-element or diagonal-by-diagonal: + +```julia +A = SparseBandedMatrix{Float64}(undef, 5, 5) +A[1,1] = 5.0 +setdiagonal!(A, [3.0, 4.0, 5.0], true) +``` + +#### Pre-specified Diagonals + +Create a matrix with known diagonals upfront: + +```julia +# Create a 6×6 matrix with two diagonals +# First diagonal at index 1 with values [3.0] +# Second diagonal at index 8 with values [-2.0, 5.0, 1.0, 3.0] +B = SparseBandedMatrix{Float64}([1, 8], [[3.0], [-2.0, 5.0, 1.0, 3.0]], 6, 6) +``` + +### Matrix Operations + +`SparseBandedMatrix` supports standard Julia array operations: + +```julia +# Matrix-vector multiplication +v = rand(5) +result = A * v + +# Matrix-matrix multiplication +B = rand(5, 3) +C = A * B + +# Element access +value = A[2, 3] +A[2, 3] = 7.0 + +# Size queries +rows, cols = size(A) +``` + +### Using `setdiagonal!` + +The `setdiagonal!` function efficiently sets an entire diagonal: + +```julia A = SparseBandedMatrix{Float64}(undef, 5, 5) -A[1,1] = 5 -setdiagonal!(A, [3,4,5], true) # sets the third diagonal from the bottom to have the values 3, 4, and 5 -B = SparseBandedMatrix{Float64}([1, 8], [[3], [-2, 5, 1, 3]], 6, 6) +# Set a lower diagonal (third from bottom) +setdiagonal!(A, [3.0, 4.0, 5.0], true) + +# Set an upper diagonal +setdiagonal!(A, [1.0, 2.0], false) ``` +## Performance Considerations + +- Use `setdiagonal!` when possible instead of setting individual elements +- Pre-specify diagonals during construction when you know them in advance +- The package is optimized for multiplication operations, making it ideal for iterative solvers + ## API Reference ```@docs diff --git a/src/SparseBandedMatrices.jl b/src/SparseBandedMatrices.jl index 3f6d89a..f772a13 100644 --- a/src/SparseBandedMatrices.jl +++ b/src/SparseBandedMatrices.jl @@ -2,6 +2,55 @@ module SparseBandedMatrices using LinearAlgebra, .Threads +""" + SparseBandedMatrix{T} <: AbstractMatrix{T} + +A sparse matrix type optimized for storing and multiplying banded matrices where only +certain diagonals contain non-zero elements. This structure is particularly efficient +for matrices arising from finite difference discretizations and butterfly factorizations. + +The matrix stores only the non-zero diagonals rather than all elements, providing +significant memory savings and computational advantages for matrix operations. + +# Constructors + + SparseBandedMatrix{T}(undef, N, M) + +Create an uninitialized `N × M` sparse banded matrix with element type `T`. +Initially contains no diagonals; diagonals are added when elements are set. + +# Arguments +- `T`: Element type (e.g., `Float64`, `ComplexF64`) +- `N`: Number of rows +- `M`: Number of columns + + SparseBandedMatrix{T}(ind_vals, diag_vals, N, M) + +Create an `N × M` sparse banded matrix with specified diagonals. + +# Arguments +- `ind_vals`: Vector of diagonal indices indicating which diagonals contain non-zero elements +- `diag_vals`: Vector of vectors containing the values for each diagonal +- `N`: Number of rows +- `M`: Number of columns + +# Examples + +```julia +# Create an empty 5×5 matrix +A = SparseBandedMatrix{Float64}(undef, 5, 5) +A[1,1] = 5.0 + +# Create a matrix with specified diagonals +B = SparseBandedMatrix{Float64}([1, 8], [[3.0], [-2.0, 5.0, 1.0, 3.0]], 6, 6) +``` + +# Implementation Notes + +Diagonals are stored internally with an indexing scheme where for an N×M matrix, +diagonal indices range from 1 to N+M-1. The storage is optimized for fast matrix +multiplication operations used in butterfly factorizations. +""" struct SparseBandedMatrix{T} <: AbstractMatrix{T} size :: Tuple{Int, Int} indices :: Vector{Int} @@ -60,8 +109,39 @@ function Base.setindex!(M :: SparseBandedMatrix{T}, val, i :: Int, j :: Int, I : val end - function setdiagonal!(M :: SparseBandedMatrix{T}, diagvals, lower :: Bool) where T - rows, cols = size(M) +""" + setdiagonal!(M::SparseBandedMatrix{T}, diagvals, lower::Bool) where T + +Set a complete diagonal of the sparse banded matrix `M` with the values in `diagvals`. + +This function efficiently sets an entire diagonal at once, which is more efficient than +setting individual elements when initializing banded matrices. + +# Arguments +- `M`: The sparse banded matrix to modify +- `diagvals`: Vector of values to assign to the diagonal. Length must not exceed the number of rows in `M` +- `lower`: If `true`, sets a lower diagonal; if `false`, sets an upper diagonal + +# Returns +Returns `diagvals` unchanged. + +# Examples + +```julia +A = SparseBandedMatrix{Float64}(undef, 5, 5) + +# Set the third diagonal from the bottom +setdiagonal!(A, [3.0, 4.0, 5.0], true) + +# Set an upper diagonal +setdiagonal!(A, [1.0, 2.0], false) +``` + +# Throws +- `ErrorException`: If `length(diagvals) > rows` (diagonal is too large for the matrix) +""" +function setdiagonal!(M :: SparseBandedMatrix{T}, diagvals, lower :: Bool) where T + rows, cols = size(M) if length(diagvals) > rows error("size of diagonal is too big for the matrix") end @@ -74,7 +154,7 @@ function Base.setindex!(M :: SparseBandedMatrix{T}, val, i :: Int, j :: Int, I : ind = searchsortedfirst(M.indices, wanted_ind) if (ind > length(M.indices) || M.indices[ind] != wanted_ind) insert!(M.indices, ind, wanted_ind) - insert!(M.diags, ind, diagvals isa Vector{T} ? diagvals : convert(Vector{T}, diagvals)::Vector{T}) + insert!(M.diags, ind, diagvals isa Vector{T} ? diagvals : convert(Vector{T}, diagvals)::Vector{T}) else for i in eachindex(diagvals) M.diags[ind][i] = diagvals[i] isa T ? diagvals[i] : convert(T, diagvals[i])::T