Improve RFP coverage and some unrelated changes to general eigenvalue problem

Author: andreasnoack

src/eigenGeneral.jl

@@ -76,8 +76,8 @@ module EigenGeneral
         iend = n
         HH = H.data
         τ = Rotation(Givens{T}[])
-        @inbounds begin
-        while true
+
+        @inbounds while true
             # Determine if the matrix splits. Find lowest positioned subdiagonal "zero"
             for istart = iend - 1:-1:1
                 # debug && @printf("istart: %6d, iend %6d\n", istart, iend)
@@ -120,8 +120,8 @@ module EigenGeneral
                     debug && @printf("Single shift! subdiagonal is: %10.3e\n", HH[iend, iend - 1])
 
                     # Calculate the Wilkinson shift
-                    λ1 = 0.5*(t + sqrt(t*t - 4d))
-                    λ2 = 0.5*(t - sqrt(t*t - 4d))
+                    λ1 = (t + sqrt(t*t - 4d))/2
+                    λ2 = (t - sqrt(t*t - 4d))/2
                     σ = abs(Hmm - λ1) < abs(Hmm - λ2) ? λ1 : λ2
 
                     # Run a bulge chase
@@ -135,10 +135,10 @@ module EigenGeneral
             end
             if iend <= 2 break end
         end
-        end
+
         return Schur{T,typeof(HH)}(HH, τ)
     end
-    schurfact!(A::StridedMatrix; tol = eps(), debug = false) = schurfact!(hessfact!(A), tol = tol, debug = debug)
+    schurfact!(A::StridedMatrix; tol = eps(float(real(eltype(A)))), debug = false) = schurfact!(hessfact!(A), tol = tol, debug = debug)
 
     function singleShiftQR!(HH::StridedMatrix, τ::Rotation, shift::Number, istart::Integer, iend::Integer)
         m = size(HH, 1)

src/lapack.jl

@@ -2,7 +2,7 @@ module LAPACK2
 
     using Base.LinAlg: BlasInt, chkstride1, LAPACKException
     using Base.LinAlg.BLAS: @blasfunc
-    using Base.LinAlg.LAPACK: chkuplo
+    using Base.LinAlg.LAPACK: chkdiag, chkuplo
 
     # LAPACK wrappers
     import Base.BLAS.@blasfunc

src/rectfullpacked.jl

@@ -114,6 +114,8 @@ function Base.size(A::TriangularRFP)
     return (n, n)
 end
 
+Base.copy(A::TriangularRFP) = TriangularRFP(copy(A.data), A.transr, A.uplo)
+
 function Base.full(A::TriangularRFP)
     C = LAPACK2.tfttr!(A.transr, A.uplo, A.data)
     if A.uplo == 'U'
@@ -123,6 +125,9 @@ function Base.full(A::TriangularRFP)
     end
 end
 
+Base.LinAlg.inv!(A::TriangularRFP) = TriangularRFP(LAPACK2.tftri!(A.transr, A.uplo, 'N', A.data), A.transr, A.uplo)
+Base.LinAlg.inv(A::TriangularRFP)  = Base.LinAlg.inv!(copy(A))
+
 struct CholeskyRFP{T<:BlasFloat} <: Factorization{T}
     data::Vector{T}
     transr::Char
@@ -139,6 +144,6 @@ Base.copy(F::CholeskyRFP{T}) where T = CholeskyRFP{T}(copy(F.data), F.transr, F.
 \(A::CholeskyRFP, B::StridedVecOrMat) = LAPACK2.pftrs!(A.transr, A.uplo, A.data, copy(B))
 \(A::HermitianRFP, B::StridedVecOrMat) = cholfact(A)\B
 
-inv!(A::CholeskyRFP) = HermitianRFP(LAPACK2.pftri!(A.transr, A.uplo, A.data), A.transr, A.uplo)
-Base.LinAlg.inv(A::CholeskyRFP)  = inv!(copy(A))
-Base.LinAlg.inv(A::HermitianRFP) = inv!(cholfact(A))
+Base.LinAlg.inv!(A::CholeskyRFP) = HermitianRFP(LAPACK2.pftri!(A.transr, A.uplo, A.data), A.transr, A.uplo)
+Base.LinAlg.inv(A::CholeskyRFP)  = Base.LinAlg.inv!(copy(A))
+Base.LinAlg.inv(A::HermitianRFP) = Base.LinAlg.inv!(cholfact(A))

test/rectfullpacked.jl

@@ -66,5 +66,6 @@ import LinearAlgebra: Ac_mul_A_RFP, TriangularRFP
 
         @test_broken A ≈ A_RFP
         @test        A ≈ full(A_RFP)
+        @test   inv(A) ≈ full(inv(A_RFP))
     end
 end