1+ # !/usr/bin/env julia
2+
3+ # Simple benchmark for SymTridiagonal LDLtFactorization improvement
4+ # This demonstrates the performance improvement from issue #732
5+
6+ using LinearSolve
7+ using LinearAlgebra
8+ using Printf
9+
10+ println (" =" ^ 80 )
11+ println (" SymTridiagonal Performance Benchmark (Issue #732)" )
12+ println (" =" ^ 80 )
13+ println ()
14+
15+ # Create test matrix from issue #732
16+ k = 10000
17+ ρ = 0.95
18+ A_tri = SymTridiagonal (ones (k) .+ ρ^ 2 , - ρ * ones (k- 1 ))
19+ b = rand (k)
20+
21+ println (" Matrix size: $k × $k (example from issue #732)" )
22+ println ()
23+
24+ # Time native \ operator
25+ println (" 1. Native \\ operator:" )
26+ t_native = @elapsed begin
27+ for _ in 1 : 10
28+ x = A_tri \ b
29+ end
30+ end
31+ println (" Average time: $(round (t_native/ 10 * 1e6 , digits= 1 )) μs" )
32+
33+ # Time LinearSolve with explicit LDLtFactorization
34+ println (" \n 2. LinearSolve with explicit LDLtFactorization:" )
35+ prob = LinearProblem (A_tri, b)
36+ t_ldlt = @elapsed begin
37+ for _ in 1 : 10
38+ x = solve (prob, LDLtFactorization ())
39+ end
40+ end
41+ println (" Average time: $(round (t_ldlt/ 10 * 1e3 , digits= 2 )) ms" )
42+
43+ # Show the factorization caching benefit
44+ println (" \n 3. Factorization caching test (100 solves with different RHS):" )
45+
46+ # Without caching (native)
47+ println (" a) Native \\ (no caching):" )
48+ t_no_cache = @elapsed begin
49+ for _ in 1 : 100
50+ b_new = rand (k)
51+ x = A_tri \ b_new
52+ end
53+ end
54+ println (" Total time: $(round (t_no_cache * 1e3 , digits= 1 )) ms" )
55+ println (" Time per solve: $(round (t_no_cache/ 100 * 1e6 , digits= 1 )) μs" )
56+
57+ # With caching (LinearSolve)
58+ println (" \n b) LinearSolve with cached LDLtFactorization:" )
59+ t_with_cache = @elapsed begin
60+ cache = init (prob, LDLtFactorization ())
61+ solve! (cache) # First solve (includes factorization)
62+ for _ in 1 : 99
63+ cache. b = rand (k)
64+ solve! (cache) # Subsequent solves (reuse factorization)
65+ end
66+ end
67+ println (" Total time: $(round (t_with_cache * 1e3 , digits= 1 )) ms" )
68+ println (" Time per solve: $(round (t_with_cache/ 100 * 1e3 , digits= 2 )) ms" )
69+
70+ speedup = t_no_cache / t_with_cache
71+ println (" \n Speedup with caching: $(round (speedup, digits= 1 )) x" )
72+
73+ # Verify the default algorithm selection
74+ println (" \n 4. Default algorithm selection:" )
75+ default_alg = LinearSolve. defaultalg (A_tri, b, LinearSolve. OperatorAssumptions (true ))
76+ if default_alg isa LinearSolve. DefaultLinearSolver
77+ if default_alg. alg == LinearSolve. DefaultAlgorithmChoice. LDLtFactorization
78+ println (" ✓ LDLtFactorization is selected by default for SymTridiagonal" )
79+ else
80+ println (" ✗ Wrong algorithm selected: $(default_alg. alg) " )
81+ end
82+ else
83+ println (" Algorithm type: $(typeof (default_alg)) " )
84+ end
85+
86+ println (" \n " * " =" ^ 80 )
87+ println (" Summary:" )
88+ println (" - LDLtFactorization works correctly with SymTridiagonal matrices" )
89+ println (" - Caching provides significant speedup for multiple solves" )
90+ println (" - Default algorithm selection has been updated to use LDLtFactorization" )
91+ println (" =" ^ 80 )
0 commit comments