Documentation improvements: Add comprehensive docstrings and examples

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. 
+## 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 

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

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