Skip to content

Commit c303ed2

Browse files
ChrisRackauckasclaude
ChrisRackauckas
and
claude
committed
Add documentation for abstract types and algorithm selection logic
- Document complete abstract type hierarchy with detailed explanations - Add comprehensive docstring for defaultalg function with algorithm selection logic - Document needs_concrete_A trait function with usage examples - Provide detailed characteristics and use cases for each abstract type - Include cross-references between related algorithm types 🤖 Generated with [Claude Code](https://claude.ai/code) Co-Authored-By: Claude <[email protected]>
1 parent 4c8d145 commit c303ed2

File tree

2 files changed

+207
-0
lines changed

2 files changed

+207
-0
lines changed

‎src/LinearSolve.jl

Lines changed: 147 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -60,15 +60,162 @@ end
6060

6161
@reexport using SciMLBase
6262

63+
"""
64+
SciMLLinearSolveAlgorithm <: SciMLBase.AbstractLinearAlgorithm
65+
66+
The root abstract type for all linear solver algorithms in LinearSolve.jl.
67+
All concrete linear solver implementations should inherit from one of the
68+
specialized subtypes rather than directly from this type.
69+
70+
This type integrates with the SciMLBase ecosystem, providing a consistent
71+
interface for linear algebra operations across the Julia scientific computing
72+
ecosystem.
73+
"""
6374
abstract type SciMLLinearSolveAlgorithm <: SciMLBase.AbstractLinearAlgorithm end
75+
76+
"""
77+
AbstractFactorization <: SciMLLinearSolveAlgorithm
78+
79+
Abstract type for linear solvers that work by computing a matrix factorization.
80+
These algorithms typically decompose the matrix `A` into a product of simpler
81+
matrices (e.g., `A = LU`, `A = QR`, `A = LDL'`) and then solve the system
82+
using forward/backward substitution.
83+
84+
## Characteristics
85+
- Requires concrete matrix representation (`needs_concrete_A() = true`)
86+
- Typically efficient for multiple solves with the same matrix
87+
- Generally provides high accuracy for well-conditioned problems
88+
- Memory requirements depend on the specific factorization type
89+
90+
## Subtypes
91+
- `AbstractDenseFactorization`: For dense matrix factorizations
92+
- `AbstractSparseFactorization`: For sparse matrix factorizations
93+
94+
## Examples of concrete subtypes
95+
- `LUFactorization`, `QRFactorization`, `CholeskyFactorization`
96+
- `UMFPACKFactorization`, `KLUFactorization`
97+
"""
6498
abstract type AbstractFactorization <: SciMLLinearSolveAlgorithm end
99+
100+
"""
101+
AbstractSparseFactorization <: AbstractFactorization
102+
103+
Abstract type for factorization-based linear solvers optimized for sparse matrices.
104+
These algorithms take advantage of sparsity patterns to reduce memory usage and
105+
computational cost compared to dense factorizations.
106+
107+
## Characteristics
108+
- Optimized for matrices with many zero entries
109+
- Often use specialized pivoting strategies to preserve sparsity
110+
- May reorder rows/columns to minimize fill-in during factorization
111+
- Typically more memory-efficient than dense methods for sparse problems
112+
113+
## Examples of concrete subtypes
114+
- `UMFPACKFactorization`: General sparse LU with partial pivoting
115+
- `KLUFactorization`: Sparse LU optimized for circuit simulation
116+
- `CHOLMODFactorization`: Sparse Cholesky for positive definite systems
117+
- `SparspakFactorization`: Envelope/profile method for sparse systems
118+
"""
65119
abstract type AbstractSparseFactorization <: AbstractFactorization end
120+
121+
"""
122+
AbstractDenseFactorization <: AbstractFactorization
123+
124+
Abstract type for factorization-based linear solvers optimized for dense matrices.
125+
These algorithms assume the matrix has no particular sparsity structure and use
126+
dense linear algebra routines (typically from BLAS/LAPACK) for optimal performance.
127+
128+
## Characteristics
129+
- Optimized for matrices with few zeros or no sparsity structure
130+
- Leverage highly optimized BLAS/LAPACK routines when available
131+
- Generally provide excellent performance for moderately-sized dense problems
132+
- Memory usage scales as O(n²) with matrix size
133+
134+
## Examples of concrete subtypes
135+
- `LUFactorization`: Dense LU with partial pivoting (via LAPACK)
136+
- `QRFactorization`: Dense QR factorization for overdetermined systems
137+
- `CholeskyFactorization`: Dense Cholesky for symmetric positive definite matrices
138+
- `BunchKaufmanFactorization`: For symmetric indefinite matrices
139+
"""
66140
abstract type AbstractDenseFactorization <: AbstractFactorization end
141+
142+
"""
143+
AbstractKrylovSubspaceMethod <: SciMLLinearSolveAlgorithm
144+
145+
Abstract type for iterative linear solvers based on Krylov subspace methods.
146+
These algorithms solve linear systems by iteratively building an approximation
147+
from a sequence of Krylov subspaces, without requiring explicit matrix factorization.
148+
149+
## Characteristics
150+
- Does not require concrete matrix representation (`needs_concrete_A() = false`)
151+
- Only needs matrix-vector products `A*v` (can work with operators/functions)
152+
- Memory usage typically O(n) or O(kn) where k << n
153+
- Convergence depends on matrix properties (condition number, eigenvalue distribution)
154+
- Often benefits significantly from preconditioning
155+
156+
## Advantages
157+
- Low memory requirements for large sparse systems
158+
- Can handle matrix-free operators (functions that compute `A*v`)
159+
- Often the only feasible approach for very large systems
160+
- Can exploit matrix structure through specialized operators
161+
162+
## Examples of concrete subtypes
163+
- `GMRESIteration`: Generalized Minimal Residual method
164+
- `CGIteration`: Conjugate Gradient (for symmetric positive definite systems)
165+
- `BiCGStabLIteration`: Bi-Conjugate Gradient Stabilized
166+
- Wrapped external iterative solvers (KrylovKit.jl, IterativeSolvers.jl)
167+
"""
67168
abstract type AbstractKrylovSubspaceMethod <: SciMLLinearSolveAlgorithm end
169+
170+
"""
171+
AbstractSolveFunction <: SciMLLinearSolveAlgorithm
172+
173+
Abstract type for linear solvers that wrap custom solving functions or
174+
provide direct interfaces to specific solve methods. These provide flexibility
175+
for integrating custom algorithms or simple solve strategies.
176+
177+
## Characteristics
178+
- Does not require concrete matrix representation (`needs_concrete_A() = false`)
179+
- Provides maximum flexibility for custom solving strategies
180+
- Can wrap external solver libraries or implement specialized algorithms
181+
- Performance and stability depend entirely on the wrapped implementation
182+
183+
## Examples of concrete subtypes
184+
- `LinearSolveFunction`: Wraps arbitrary user-defined solve functions
185+
- `DirectLdiv!`: Direct application of the `\\` operator
186+
"""
68187
abstract type AbstractSolveFunction <: SciMLLinearSolveAlgorithm end
69188

70189
# Traits
71190

191+
"""
192+
needs_concrete_A(alg) -> Bool
193+
194+
Trait function that determines whether a linear solver algorithm requires
195+
a concrete matrix representation or can work with abstract operators.
196+
197+
## Arguments
198+
- `alg`: A linear solver algorithm instance
199+
200+
## Returns
201+
- `true`: Algorithm requires a concrete matrix (e.g., for factorization)
202+
- `false`: Algorithm can work with abstract operators (e.g., matrix-free methods)
203+
204+
## Usage
205+
This trait is used internally by LinearSolve.jl to optimize algorithm dispatch
206+
and determine when matrix operators need to be converted to concrete arrays.
207+
208+
## Algorithm-Specific Behavior
209+
- `AbstractFactorization`: `true` (needs explicit matrix entries for factorization)
210+
- `AbstractKrylovSubspaceMethod`: `false` (only needs matrix-vector products)
211+
- `AbstractSolveFunction`: `false` (depends on the wrapped function's requirements)
212+
213+
## Example
214+
```julia
215+
needs_concrete_A(LUFactorization()) # true
216+
needs_concrete_A(GMRESIteration()) # false
217+
```
218+
"""
72219
needs_concrete_A(alg::AbstractFactorization) = true
73220
needs_concrete_A(alg::AbstractSparseFactorization) = true
74221
needs_concrete_A(alg::AbstractKrylovSubspaceMethod) = false

‎src/default.jl

Lines changed: 60 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -47,6 +47,66 @@ function __setfield!(cache::DefaultLinearSolverInit,
4747
setfield!(cache, :QRFactorizationPivoted, v)
4848
end
4949

50+
"""
51+
defaultalg(A, b, assumptions::OperatorAssumptions)
52+
53+
Select the most appropriate linear solver algorithm based on the matrix `A`,
54+
right-hand side `b`, and operator assumptions. This is the core algorithm
55+
selection logic used by LinearSolve.jl's automatic algorithm choice.
56+
57+
## Arguments
58+
- `A`: The matrix operator (can be a matrix, factorization, or abstract operator)
59+
- `b`: The right-hand side vector
60+
- `assumptions`: Operator assumptions including square matrix flag and conditioning
61+
62+
## Returns
63+
A `DefaultLinearSolver` instance configured with the most appropriate algorithm choice,
64+
or a specific algorithm instance for certain special cases.
65+
66+
## Algorithm Selection Logic
67+
68+
The function uses a hierarchy of dispatch rules based on:
69+
70+
1. **Matrix Type**: Special handling for structured matrices (Diagonal, Tridiagonal, etc.)
71+
2. **Matrix Properties**: Square vs. rectangular, sparse vs. dense
72+
3. **Hardware**: GPU vs. CPU arrays
73+
4. **Conditioning**: Well-conditioned vs. ill-conditioned systems
74+
5. **Size**: Small vs. large matrices for performance optimization
75+
76+
## Common Algorithm Choices
77+
78+
- **Diagonal matrices**: `DiagonalFactorization` for optimal O(n) performance
79+
- **Tridiagonal/Bidiagonal**: Direct methods or specialized factorizations
80+
- **Dense matrices**: LU, QR, or Cholesky based on structure and conditioning
81+
- **Sparse matrices**: Specialized sparse factorizations (UMFPACK, KLU, etc.)
82+
- **GPU arrays**: QR or LU factorizations optimized for GPU computation
83+
- **Abstract operators**: Krylov methods (GMRES, CRAIGMR, LSMR)
84+
- **Symmetric positive definite**: Cholesky factorization
85+
- **Symmetric indefinite**: Bunch-Kaufman factorization
86+
87+
## Examples
88+
89+
```julia
90+
# Dense square matrix - typically chooses LU
91+
A = rand(100, 100)
92+
b = rand(100)
93+
alg = defaultalg(A, b, OperatorAssumptions(true))
94+
95+
# Overdetermined system - typically chooses QR
96+
A = rand(100, 50)
97+
b = rand(100)
98+
alg = defaultalg(A, b, OperatorAssumptions(false))
99+
100+
# Diagonal matrix - chooses diagonal factorization
101+
A = Diagonal(rand(100))
102+
alg = defaultalg(A, b, OperatorAssumptions(true))
103+
```
104+
105+
## Notes
106+
This function is primarily used internally by `solve(::LinearProblem)` when no
107+
explicit algorithm is provided. For manual algorithm selection, users can
108+
directly instantiate specific algorithm types.
109+
"""
50110
# Legacy fallback
51111
# For SciML algorithms already using `defaultalg`, all assume square matrix.
52112
defaultalg(A, b) = defaultalg(A, b, OperatorAssumptions(true))

0 commit comments

Comments
 (0)