Fix slow triangular matrix solves by using DirectLdiv\! (#672)

* Fix slow triangular matrix solves by using DirectLdiv\! directly

Addresses issue #671 by changing the default algorithm selection for
triangular matrices to use DirectLdiv\!() directly instead of routing
through DefaultLinearSolver, which has expensive initialization overhead.

Changes:
- SymTridiagonal: LDLtFactorization → DirectLdiv\!()
- Tridiagonal: LUFactorization → DirectLdiv\!()
- Bidiagonal: DirectLdiv\!() (consistent with others)

This approach uses type-inferred dispatch to avoid the DefaultLinearSolver's
overhead of initializing cache for all ~20 algorithms when only one is needed.

Performance improvement: ~70x faster for large triangular matrices,
now matching or exceeding native \ operator performance.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <[email protected]>

* Update src/default.jl

* Update src/default.jl

* Update src/default.jl

---------

Co-authored-by: ChrisRackauckas <[email protected]>
Co-authored-by: Claude <[email protected]>

Author: 3 people

src/default.jl

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