Use LDLtFactorization for SymTridiagonal matrices (fixes #732) (#735)

* Use LDLtFactorization for SymTridiagonal matrices

- Changed default algorithm selection for SymTridiagonal to use LDLtFactorization
- This caches the LDLT factorization for reuse in subsequent solves
- Addresses performance issue where ldlt was being recomputed on every solve

* Add comprehensive tests for SymTridiagonal with LDLtFactorization

- Test explicit LDLtFactorization with SymTridiagonal matrices
- Verify default algorithm selects LDLtFactorization for SymTridiagonal
- Test that factorization is cached and reused across solves
- Include performance test case from issue #732

* Add performance benchmark demonstrating the fix

- Shows LDLtFactorization working correctly with SymTridiagonal
- Demonstrates 5.8x speedup with factorization caching
- Verifies default algorithm selection

* Delete benchmark_simple.jl

* Update test/basictests.jl

---------

Co-authored-by: Claude <noreply@anthropic.com>
Co-authored-by: Christopher Rackauckas <accounts@chrisrackauckas.com>

Author: 3 people

src/default.jl

@@ -145,11 +145,7 @@ function defaultalg(A::Tridiagonal, b, assump::OperatorAssumptions{Bool})
 end
 
 function defaultalg(A::SymTridiagonal, b, ::OperatorAssumptions{Bool})
-    @static if VERSION>=v"1.11"
-          DirectLdiv!()
-      else
-          DefaultLinearSolver(DefaultAlgorithmChoice.LUFactorization)
-      end
+    DefaultLinearSolver(DefaultAlgorithmChoice.LDLtFactorization)
 end
 function defaultalg(A::Bidiagonal, b, ::OperatorAssumptions{Bool})
     @static if VERSION>=v"1.11"

test/basictests.jl

@@ -275,6 +275,37 @@ end
         # ws sizes.
     end
 
+    @testset "SymTridiagonal with LDLtFactorization" begin
+        # Test that LDLtFactorization works correctly with SymTridiagonal
+        # and that the default algorithm correctly selects it
+        k = 100
+        ρ = 0.95
+        A_tri = SymTridiagonal(ones(k) .+ ρ^2, -ρ * ones(k-1))
+        b = rand(k)
+        
+        # Test with explicit LDLtFactorization
+        prob_tri = LinearProblem(A_tri, b)
+        sol = solve(prob_tri, LDLtFactorization())
+        @test A_tri * sol.u ≈ b
+        
+        # Test that default algorithm uses LDLtFactorization for SymTridiagonal
+        default_alg = LinearSolve.defaultalg(A_tri, b, OperatorAssumptions(true))
+        @test default_alg isa LinearSolve.DefaultLinearSolver
+        @test default_alg.alg == LinearSolve.DefaultAlgorithmChoice.LDLtFactorization
+        
+        # Test that the factorization is cached and reused
+        cache = init(prob_tri, LDLtFactorization())
+        sol1 = solve!(cache)
+        @test A_tri * sol1.u ≈ b
+        @test !cache.isfresh  # Cache should not be fresh after first solve
+        
+        # Solve again with same matrix to ensure cache is reused
+        cache.b = rand(k)  # Change RHS
+        sol2 = solve!(cache)
+        @test A_tri * sol2.u ≈ cache.b
+        @test !cache.isfresh  # Cache should still not be fresh
+    end
+
     test_algs = [
         LUFactorization(),
         QRFactorization(),