Use BLAS.trsm! instead of LAPACK.trtrs! in left-triangular solves (#1194)

Co-authored-by: Alexis Montoison <[email protected]>

Author: dkarrasch

Diff for: src/blas.jl

@@ -84,7 +84,8 @@ export
     trsm!,
     trsm
 
-using ..LinearAlgebra: libblastrampoline, BlasReal, BlasComplex, BlasFloat, BlasInt, DimensionMismatch, checksquare, chkstride1
+using ..LinearAlgebra: libblastrampoline, BlasReal, BlasComplex, BlasFloat, BlasInt,
+    DimensionMismatch, checksquare, chkstride1, SingularException
 
 include("lbt.jl")
 
@@ -1369,6 +1370,11 @@ for (fname, elty) in ((:dtrsv_,:Float64),
                 throw(DimensionMismatch(lazy"size of A is $n != length(x) = $(length(x))"))
             end
             chkstride1(A)
+            if diag == 'N'
+                for i in 1:n
+                    iszero(A[i,i]) && throw(SingularException(i))
+                end
+            end
             px, stx = vec_pointer_stride(x, ArgumentError("input vector with 0 stride is not allowed"))
             GC.@preserve x ccall((@blasfunc($fname), libblastrampoline), Cvoid,
                 (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt},
@@ -2217,6 +2223,11 @@ for (mmname, smname, elty) in
             end
             chkstride1(A)
             chkstride1(B)
+            if diag == 'N'
+                for i in 1:k
+                    iszero(A[i,i]) && throw(SingularException(i))
+                end
+            end
             ccall((@blasfunc($smname), libblastrampoline), Cvoid,
                    (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{UInt8},
                     Ref{BlasInt}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty},

Diff for: src/triangular.jl

@@ -1223,11 +1223,13 @@ function generic_mattrimul!(C::StridedMatrix{T}, uploc, isunitc, tfun::Function,
     end
 end
 # division
-function generic_trimatdiv!(C::StridedVecOrMat{T}, uploc, isunitc, tfun::Function, A::StridedMatrix{T}, B::AbstractVecOrMat{T}) where {T<:BlasFloat}
+generic_trimatdiv!(C::StridedVector{T}, uploc, isunitc, tfun::Function, A::StridedMatrix{T}, B::AbstractVector{T}) where {T<:BlasFloat} =
+    BLAS.trsv!(uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, A, C === B ? C : copyto!(C, B))
+function generic_trimatdiv!(C::StridedMatrix{T}, uploc, isunitc, tfun::Function, A::StridedMatrix{T}, B::AbstractMatrix{T}) where {T<:BlasFloat}
     if stride(C,1) == stride(A,1) == 1
-        LAPACK.trtrs!(uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, A, C === B ? C : copyto!(C, B))
+        BLAS.trsm!('L', uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, one(T), A, C === B ? C : copyto!(C, B))
     else # incompatible with LAPACK
-        @invoke generic_trimatdiv!(C::AbstractVecOrMat, uploc, isunitc, tfun::Function, A::AbstractMatrix, B::AbstractVecOrMat)
+        @invoke generic_trimatdiv!(C::AbstractVecOrMat, uploc, isunitc, tfun::Function, A::AbstractMatrix, B::AbstractMatrix)
     end
 end
 function generic_mattridiv!(C::StridedMatrix{T}, uploc, isunitc, tfun::Function, A::AbstractMatrix{T}, B::StridedMatrix{T}) where {T<:BlasFloat}

Diff for: test/testtriag.jl

@@ -493,6 +493,8 @@ function test_triangular(elty1_types)
                 @test_throws DimensionMismatch transpose(Ann) \ bm
                 if t1 == UpperTriangular || t1 == LowerTriangular
                     @test_throws SingularException ldiv!(t1(zeros(elty1, n, n)), fill(eltyB(1), n))
+                    @test_throws SingularException ldiv!(t1(zeros(elty1, n, n)), fill(eltyB(1), n, 2))
+                    @test_throws SingularException rdiv!(fill(eltyB(1), n, n), t1(zeros(elty1, n, n)))
                 end
                 @test B / A1 ≈ B / M1
                 @test B / transpose(A1) ≈ B / transpose(M1)

Diff for: test/triangular.jl

@@ -886,8 +886,13 @@ end
     end
 end
 
-@testset "(l/r)mul! and (l/r)div! for non-contiguous matrices" begin
+@testset "(l/r)mul! and (l/r)div! for non-contiguous arrays" begin
     U = UpperTriangular(reshape(collect(3:27.0),5,5))
+    b = float.(1:10)
+    b2 = copy(b); b2v = view(b2, 1:2:9); b2vc = copy(b2v)
+    @test lmul!(U, b2v) == lmul!(U, b2vc)
+    b2 = copy(b); b2v = view(b2, 1:2:9); b2vc = copy(b2v)
+    @test ldiv!(U, b2v) ≈ ldiv!(U, b2vc)
     B = float.(collect(reshape(1:100, 10,10)))
     B2 = copy(B); B2v = view(B2, 1:2:9, 1:5); B2vc = copy(B2v)
     @test lmul!(U, B2v) == lmul!(U, B2vc)