fix standard real form logic and make it the default again

src/GenericSchur.jl

@@ -256,16 +256,23 @@ function singleShiftQR!(HH::StridedMatrix{T}, Z, shift::Number, istart::Integer,
     return HH
 end
 
+const _STANDARDIZE_DEFAULT = true
+
 function _gschur!(H::HessenbergFactorization{T}, Z=nothing;
                   tol = eps(real(T)), debug = false, shiftmethod = :Francis,
-                  maxiter = 100*size(H, 1), standardize = false,
+                  maxiter = 100*size(H, 1), standardize = _STANDARDIZE_DEFAULT,
                   kwargs...) where {T <: Real}
     n = size(H, 1)
     istart = 1
     iend = n
     HH = H.data
     τ = Rotation(Givens{T}[])
-    w = Complex{T}[]
+    w = Vector{Complex{T}}(undef, n)
+    iwcur = n
+    function putw!(w,z)
+        w[iwcur] = z
+        iwcur -= 1
+    end
 
     # iteration count
     i = 0
@@ -302,7 +309,7 @@ function _gschur!(H::HessenbergFactorization{T}, Z=nothing;
         # if block size is one we deflate
         if istart >= iend
             debug && @printf("Bottom deflation! Block size is one. New iend is %6d\n", iend - 1)
-            standardize && push!(w,HH[iend,iend])
+            standardize && putw!(w,HH[iend,iend])
             iend -= 1
 
         # and the same for a 2x2 block
@@ -313,8 +320,8 @@ function _gschur!(H::HessenbergFactorization{T}, Z=nothing;
                 debug && println("std. iend = $iend")
                 H2 = HH[iend-1:iend,iend-1:iend]
                 G2,w1,w2 = _gs2x2!(H2,iend)
-                push!(w,w1)
-                push!(w,w2)
+                putw!(w,w2)
+                putw!(w,w1)
                 lmul!(G2,view(HH,:,istart:n))
                 rmul!(view(HH,1:iend,:),G2')
                 HH[iend-1:iend,iend-1:iend] .= H2 # clean
@@ -388,15 +395,16 @@ function _gschur!(H::HessenbergFactorization{T}, Z=nothing;
         if iend <= 2
             if standardize
                 if iend == 1
-                    push!(w,HH[iend,iend])
+                    putw!(w,HH[iend,iend])
                 elseif iend == 2
-                    H2 = HH[iend-1:iend,iend-1:iend]
-                    G2,w1,w2 = _gs2x2!(H2,iend)
-                    push!(w,w1)
-                    push!(w,w2)
+                    debug && println("final std. iend = $iend")
+                    H2 = HH[1:2,1:2]
+                    G2,w1,w2 = _gs2x2!(H2,2)
+                    putw!(w,w2)
+                    putw!(w,w1)
                     lmul!(G2,HH)
                     rmul!(view(HH,1:2,:),G2')
-                    HH[iend-1:iend,iend-1:iend] .= H2 # clean
+                    HH[1:2,1:2] .= H2 # clean
                     Z === nothing || rmul!(Z,G2')
                 end
             end
@@ -479,6 +487,7 @@ function _gs2x2!(H2::StridedMatrix{T},jj) where {T <: Real}
                         τ = one(T) / sqrt(abs(b+c))
                         a = midad + p
                         d = midad - p
+                        b -= c
                         c = 0
                         cs1 = sab*τ
                         sn1 = sac*τ

test/real.jl

@@ -1,10 +1,12 @@
-const _standard = Ref(false)
+# This is the value of standardize used for the vast majority of tests:
+const _standard = Ref(true)
 
 if !_standard[]
     @info "Suppressing 2x2 block checks pending standardization fixes."
 end
 
-# this requires that tiny values have been cleaned out
+# Check that 2x2 diagonal blocks have standard form.
+# This requires that tiny values have been cleaned out.
 function checkblocks(A::StridedMatrix; debug=false)
     n = size(A,1)
     isok = true
@@ -26,12 +28,43 @@ function checkblocks(A::StridedMatrix; debug=false)
     isok
 end
 
+# This only applies to standard forms w/ co-sorted eigenvals
+function checkeigvals(A::StridedMatrix{T}, w::StridedVector, tol; debug=false
+                      ) where T <: Real
+    n = size(A,1)
+    ulp = eps(T)
+    tiny = GenericSchur.safemin(T)
+    isok = true
+    for j=1:n
+        isok &= (A[j,j] == real(w[j]))
+    end
+    if n>1
+        if A[2,1] == 0
+            isok &= (imag(w[1]) == 0)
+        end
+        if A[n,n-1] == 0
+            isok &= (imag(w[n]) == 0)
+        end
+    end
+    for j=1:n-1
+        if A[j+1,j] != 0
+            t = sqrt(abs(A[j+1,j]))*sqrt(abs(A[j,j+1]))
+            cmp = max(ulp*t,tiny)
+            isok &= (abs(imag(w[j]) - t) / cmp < tol)
+            isok &= (abs(imag(w[j+1]) + t) / cmp < tol)
+        elseif (j>1) && (A[j+1,j] == 0) && (A[j,j-1] == 0)
+            isok &= (imag(w[j]) == 0)
+        end
+    end
+    isok
+end
+
 function schurtest(A::Matrix{T}, tol; normal=false, standard=_standard[],
                    baddec=false,
                    ) where {T<:Real}
     n = size(A,1)
     ulp = eps(T)
-    if T <: BlasFloat
+    if T <: BlasFloat || (standard != GenericSchur._STANDARDIZE_DEFAULT)
         S = GenericSchur.gschur(A, standardize=standard)
     else
         # FIXME: no keywords allowed here; thanks, LinearAlgebra
@@ -52,15 +85,9 @@ function schurtest(A::Matrix{T}, tol; normal=false, standard=_standard[],
     # test 3: S.Z is orthogonal: norm(I - S.Z * S.Z') / (n * ulp) < tol
     @test norm(I - S.Z * S.Z') / (n * ulp) < tol
     # test 4: S.values are e.v. of T
-    # I thought this would be robust enough (hey, IWFM),
-    # but Travis runs on hardware that says otherwise.
-    if false # (T <: BlasFloat) && !standard
-        vLA = csort(eigvals(S.T))
-        vGS = csort(S.values)
-        rte = sqrt(eps(T))
-        @test isapprox(vGS, vLA, atol=rte, rtol=rte)
+    if standard
+        @test checkeigvals(S.T, S.values, tol)
     end
-    # TODO: in standard case test 4 should be part of checkblocks
 
     # It is tempting to check eigenvalues against LAPACK eigvals(A)
     # for suitable types, but the comparison is misleading without
@@ -250,8 +277,9 @@ for T in [Float64, Float32]
                     (verbosity[] > 1) &&
                         println("latmr: n=$n $anorm, $imode, $rcond")
                     latmr!(A,anorm,imode,rcond)
-                    schurtest(A,tols[itype],
-                              baddec=(_standard[] && (imode==6) && (j==3) && (n==32)))
+                    schurtest(A,tols[itype])
+                    # throw this in for coverage
+                    schurtest(A,tols[itype],standard=!_standard[])
                 end
             end
         end