Merge pull request #117 from JuliaLinearAlgebra/an/morecoverage

Fix and test pteqr

Author: andreasnoack

src/lapack.jl

@@ -425,31 +425,39 @@ function lahqr!(H::StridedMatrix{Float64})
 end
 
 ## Cholesky + singular values
-function spteqr!(
+function pteqr!(
     compz::Char,
     d::StridedVector{Float64},
     e::StridedVector{Float64},
     Z::StridedMatrix{Float64},
-    work::StridedVector{Float64} = Vector{Float64}(undef, 4length(d)),
+    work::StridedVector{Float64},
 )
-
     n = length(d)
     ldz = stride(Z, 2)
 
     # Checks
-    length(e) >= n - 1 || throw(DimensionMismatch("subdiagonal vector is too short"))
+    if length(e) < n - 1
+        throw(DimensionMismatch("subdiagonal vector is too short"))
+    end
+    chkstride1(d)
+    chkstride1(e)
+    chkstride1(Z)
     if compz == 'N'
-    elseif compz == in('V', 'I')
-        size(Z, 1) >= n || throw(DimensionMismatch("Z does not have enough rows"))
-        size(Z, 2) >= ldz || throw(DimensionMismatch("Z does not have enough columns"))
+    elseif compz ∈ ('V', 'I')
+        if size(Z, 1) < n
+            throw(DimensionMismatch("Z does not have enough rows"))
+        end
+        if size(Z, 2) < ldz
+            throw(DimensionMismatch("Z does not have enough columns"))
+        end
     else
         throw(ArgumentError("compz must be either 'N', 'V', or 'I'"))
     end
 
-    info = Vector{BlasInt}(undef, 1)
+    info = Ref{BlasInt}(1)
 
     ccall(
-        (@blasfunc(:dpteqr_), liblapack_name),
+        (@blasfunc(dpteqr_), liblapack_name),
         Cvoid,
         (
             Ref{UInt8},
@@ -459,7 +467,7 @@ function spteqr!(
             Ptr{Float64},
             Ref{BlasInt},
             Ptr{Float64},
-            Ptr{BlasInt},
+            Ref{BlasInt},
         ),
         compz,
         n,
@@ -471,9 +479,23 @@ function spteqr!(
         info,
     )
 
-    info[1] == 0 || throw(LAPACKException(info[1]))
+    if info[] != 0
+        throw(LAPACKException(info[]))
+    end
+
+    return d, Z
+end
 
-    d, Z
+function pteqr!(compz::Char, d::StridedVector{Float64}, e::StridedVector{Float64})
+    n = length(d)
+
+    Z = if compz == 'I'
+        Matrix{Float64}(undef, n, n)
+    else
+        Matrix{Float64}(undef, 0, 0)
+    end
+    work = Vector{Float64}(undef, 4 * n)
+    return pteqr!(compz, d, e, Z, work)
 end
 
 # Gu's dnc eigensolver

test/lapack.jl

@@ -91,7 +91,7 @@ using GenericLinearAlgebra.LAPACK2
                                                                      (Float32, Float64),
         side in ('L', 'R', 'B'),
         howmny in ('A', 'S')
- #='B', =#
+
         select = ones(Int, n)
         S, P = triu(randn(eltype, n, n)), triu(randn(eltype, n, n))
         VL, VR, m = LAPACK2.tgevc!(side, howmny, select, copy(S), copy(P))
@@ -106,4 +106,14 @@ using GenericLinearAlgebra.LAPACK2
                 sqrt(eps(eltype))
         end
     end
+
+    @testset "pteqr" begin
+        d = fill(10.0, n)
+        e = fill(1.0, n - 1)
+        vals, vecs = LAPACK2.pteqr!('I', copy(d), copy(e))
+        @test SymTridiagonal(d, e) ≈ vecs * Diagonal(vals) * vecs'
+
+        vals2, _ = LAPACK2.pteqr!('N', copy(d), copy(e))
+        @test vals ≈ vals2
+    end
 end

test/svd.jl

@@ -239,6 +239,6 @@ using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions, DoubleFloats
     @testset "Generic HessenbergQ multiplication" begin
         A = big.(randn(10, 10))
         BF = GenericLinearAlgebra.bidiagonalize!(copy(A))
-        @test (BF.rightQ'*Matrix(I, size(A)...))*BF.rightQ ≈ I
+        @test (BF.rightQ' * Matrix(I, size(A)...)) * BF.rightQ ≈ I
     end
 end