Complete native binary triangular solve for MKL and BLIS (#667)

* Complete MKL triangular solve with native LAPACK calls

Replace Julia ldiv\! fallback with direct MKL getrs\! calls for the triangular
solve portion of MKLLUFactorization. This ensures the entire LU solve process
uses native MKL LAPACK routines instead of falling back to libblastrampoline.

Changes:
- Use existing getrs\! functions that were already implemented but unused
- Handle both square and overdetermined systems with proper dimension checks
- Add proper error handling for failed factorizations
- Maintain compatibility with existing LinearCache interface

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

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

* Complete BLIS triangular solve with native LAPACK calls

Replace Julia ldiv\! fallback with direct LAPACK getrs\! calls via BLIS for the
triangular solve portion of BLISLUFactorization. This ensures the entire LU
solve process uses native LAPACK routines through BLIS instead of falling back
to libblastrampoline.

Changes:
- Use existing getrs\! functions that were already implemented but unused
- Handle both square and overdetermined systems with proper dimension checks
- Add proper error handling for failed factorizations with ReturnCode
- Add missing ReturnCode import from SciMLBase
- Maintain compatibility with existing LinearCache interface

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

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

* Update Tests.yml

---------

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

Author: 3 people

.github/workflows/Tests.yml

@@ -36,6 +36,10 @@ jobs:
           - "LinearSolveHYPRE"
           - "LinearSolvePardiso"
           - "NoPre"
+        os:
+          - ubuntu-latest
+          - macos-latest
+          - windows-latest
     uses: "SciML/.github/.github/workflows/tests.yml@v1"
     with:
       group: "${{ matrix.group }}"

ext/LinearSolveBLISExt.jl

@@ -10,6 +10,7 @@ using LinearAlgebra: BlasInt, LU
 using LinearAlgebra.LAPACK: require_one_based_indexing, chkfinite, chkstride1, 
                             @blasfunc, chkargsok
 using LinearSolve: ArrayInterface, BLISLUFactorization, @get_cacheval, LinearCache, SciMLBase
+using SciMLBase: ReturnCode
 
 const global libblis = blis_jll.blis
 const global liblapack = LAPACK_jll.liblapack
@@ -224,27 +225,27 @@ function SciMLBase.solve!(cache::LinearCache, alg::BLISLUFactorization;
         res = getrf!(A; ipiv = cacheval[1].ipiv, info = cacheval[2])
         fact = LU(res[1:3]...), res[4]
         cache.cacheval = fact
+
+        if !LinearAlgebra.issuccess(fact[1])
+            return SciMLBase.build_linear_solution(
+                alg, cache.u, nothing, cache; retcode = ReturnCode.Failure)
+        end
         cache.isfresh = false
     end
 
-    y = ldiv!(cache.u, @get_cacheval(cache, :BLISLUFactorization)[1], cache.b)
-    SciMLBase.build_linear_solution(alg, y, nothing, cache)
-
-    #=
     A, info = @get_cacheval(cache, :BLISLUFactorization)
-    LinearAlgebra.require_one_based_indexing(cache.u, cache.b)
+    require_one_based_indexing(cache.u, cache.b)
     m, n = size(A, 1), size(A, 2)
     if m > n
         Bc = copy(cache.b)
         getrs!('N', A.factors, A.ipiv, Bc; info)
-        return copyto!(cache.u, 1, Bc, 1, n)
+        copyto!(cache.u, 1, Bc, 1, n)
     else
         copyto!(cache.u, cache.b)
         getrs!('N', A.factors, A.ipiv, cache.u; info)
     end
 
-    SciMLBase.build_linear_solution(alg, cache.u, nothing, cache)
-    =#
+    SciMLBase.build_linear_solution(alg, cache.u, nothing, cache; retcode = ReturnCode.Success)
 end
 
 end

src/mkl.jl

@@ -227,22 +227,17 @@ function SciMLBase.solve!(cache::LinearCache, alg::MKLLUFactorization;
         cache.isfresh = false
     end
 
-    y = ldiv!(cache.u, @get_cacheval(cache, :MKLLUFactorization)[1], cache.b)
-    SciMLBase.build_linear_solution(alg, y, nothing, cache; retcode = ReturnCode.Success)
-
-    #=
     A, info = @get_cacheval(cache, :MKLLUFactorization)
-    LinearAlgebra.require_one_based_indexing(cache.u, cache.b)
+    require_one_based_indexing(cache.u, cache.b)
     m, n = size(A, 1), size(A, 2)
     if m > n
         Bc = copy(cache.b)
         getrs!('N', A.factors, A.ipiv, Bc; info)
-        return copyto!(cache.u, 1, Bc, 1, n)
+        copyto!(cache.u, 1, Bc, 1, n)
     else
         copyto!(cache.u, cache.b)
         getrs!('N', A.factors, A.ipiv, cache.u; info)
     end
 
-    SciMLBase.build_linear_solution(alg, cache.u, nothing, cache)
-    =#
+    SciMLBase.build_linear_solution(alg, cache.u, nothing, cache; retcode = ReturnCode.Success)
 end