From 6679a390e376f54a7e12bc0cc8fc95746e3c9ad0 Mon Sep 17 00:00:00 2001
From: Simeon Schaub <simeondavidschaub99@gmail.com>
Date: Mon, 14 Jun 2021 21:40:14 +0200
Subject: [PATCH 1/3] fix equality of QRCompactWY

Equality for `QRCompactWY` did not ignore the subdiagonal entries of
`T` leading to nondeterministic behavior. Perhaps `T` should be
directly stored as `UpperTriangular` in `QRCompactWY`, but that seems
potentially breaking.

This is pulled out from #41228, since this change should be less
controversial than the other changes there and this particular bug just
came up in ChainRules again.
---
 stdlib/LinearAlgebra/src/qr.jl             | 10 +++++
 stdlib/LinearAlgebra/test/factorization.jl | 45 ++++++++++++++++++++++
 stdlib/LinearAlgebra/test/testgroups       |  1 +
 3 files changed, 56 insertions(+)
 create mode 100644 stdlib/LinearAlgebra/test/factorization.jl

diff --git a/stdlib/LinearAlgebra/src/qr.jl b/stdlib/LinearAlgebra/src/qr.jl
index 3bf5877192410..2be6b6ff573b4 100644
--- a/stdlib/LinearAlgebra/src/qr.jl
+++ b/stdlib/LinearAlgebra/src/qr.jl
@@ -127,6 +127,16 @@ Base.iterate(S::QRCompactWY) = (S.Q, Val(:R))
 Base.iterate(S::QRCompactWY, ::Val{:R}) = (S.R, Val(:done))
 Base.iterate(S::QRCompactWY, ::Val{:done}) = nothing
 
+function Base.hash(F::QRCompactWY, h::UInt)
+    return hash(F.factors, hash(UpperTriangular(F.T), hash(QRCompactWY, h)))
+end
+function Base.:(==)(A::QRCompactWY, B::QRCompactWY)
+    return A.factors == B.factors && UpperTriangular(A.T) == UpperTriangular(B.T)
+end
+function Base.isequal(A::QRCompactWY, B::QRCompactWY)
+    return isequal(A.factors, B.factors) && isequal(UpperTriangular(A.T), UpperTriangular(B.T))
+end
+
 """
     QRPivoted <: Factorization
 
diff --git a/stdlib/LinearAlgebra/test/factorization.jl b/stdlib/LinearAlgebra/test/factorization.jl
new file mode 100644
index 0000000000000..a67871af4cdb0
--- /dev/null
+++ b/stdlib/LinearAlgebra/test/factorization.jl
@@ -0,0 +1,45 @@
+# This file is a part of Julia. License is MIT: https://julialang.org/license
+
+module TestFactorization
+using Test, LinearAlgebra
+
+@testset "equality for factorizations - $f" for f in Any[
+    bunchkaufman,
+    cholesky,
+    x -> cholesky(x, Val(true)),
+    eigen,
+    hessenberg,
+    lq,
+    lu,
+    qr,
+    x -> qr(x, ColumnNorm()),
+    svd,
+    schur,
+]
+    A = randn(3, 3)
+    A = A * A' # ensure A is pos. def. and symmetric
+    F, G = f(A), f(A)
+
+    @test F == G
+    @test isequal(F, G)
+    @test hash(F) == hash(G)
+
+    f === hessenberg && continue
+
+    # change all arrays in F to have eltype Float32
+    F = typeof(F).name.wrapper(Base.mapany(1:nfields(F)) do i
+        x = getfield(F, i)
+        return x isa AbstractArray{Float64} ? Float32.(x) : x
+    end...)
+    # round all arrays in G to the nearest Float64 representable as Float32
+    G = typeof(G).name.wrapper(Base.mapany(1:nfields(G)) do i
+        x = getfield(G, i)
+        return x isa AbstractArray{Float64} ? Float64.(Float32.(x)) : x
+    end...)
+
+    @test F == G broken=!(f === eigen || f === qr)
+    @test isequal(F, G) broken=!(f === eigen || f === qr)
+    @test hash(F) == hash(G)
+end
+
+end
diff --git a/stdlib/LinearAlgebra/test/testgroups b/stdlib/LinearAlgebra/test/testgroups
index b33dfecaa82ee..de082d8e7dce0 100644
--- a/stdlib/LinearAlgebra/test/testgroups
+++ b/stdlib/LinearAlgebra/test/testgroups
@@ -25,3 +25,4 @@ givens
 structuredbroadcast
 addmul
 ldlt
+factorization

From 72e3b410a61a48410fc7a2c4deb6a1594d1fb05e Mon Sep 17 00:00:00 2001
From: Simeon Schaub <simeondavidschaub99@gmail.com>
Date: Sun, 27 Jun 2021 23:38:16 +0200
Subject: [PATCH 2/3] fix for matrix sizes >36

---
 stdlib/LinearAlgebra/src/LinearAlgebra.jl  |  3 ++-
 stdlib/LinearAlgebra/src/qr.jl             | 28 +++++++++++++++++++---
 stdlib/LinearAlgebra/test/factorization.jl | 15 ++++++++++++
 3 files changed, 42 insertions(+), 4 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl
index ebdb78a0c63c2..5f91b4832069c 100644
--- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl
+++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl
@@ -16,7 +16,8 @@ import Base: USE_BLAS64, abs, acos, acosh, acot, acoth, acsc, acsch, adjoint, as
     setindex!, show, similar, sin, sincos, sinh, size, sqrt,
     strides, stride, tan, tanh, transpose, trunc, typed_hcat, vec
 using Base: IndexLinear, promote_op, promote_typeof,
-    @propagate_inbounds, @pure, reduce, typed_vcat, require_one_based_indexing
+    @propagate_inbounds, @pure, reduce, typed_vcat, require_one_based_indexing,
+    splat
 using Base.Broadcast: Broadcasted, broadcasted
 import Libdl
 
diff --git a/stdlib/LinearAlgebra/src/qr.jl b/stdlib/LinearAlgebra/src/qr.jl
index 2be6b6ff573b4..1ebf0231fb1ad 100644
--- a/stdlib/LinearAlgebra/src/qr.jl
+++ b/stdlib/LinearAlgebra/src/qr.jl
@@ -127,14 +127,36 @@ Base.iterate(S::QRCompactWY) = (S.Q, Val(:R))
 Base.iterate(S::QRCompactWY, ::Val{:R}) = (S.R, Val(:done))
 Base.iterate(S::QRCompactWY, ::Val{:done}) = nothing
 
+# returns upper triangular views of all non-undef values of `qr(A).T`:
+#
+# julia> sparse(qr(A).T .== qr(A).T)
+# 36×100 SparseMatrixCSC{Bool, Int64} with 1767 stored entries:
+# ⠙⢿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⠙⢿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⠙⢿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
+# ⠀⠀⠙⢿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⠀⠀⠙⢿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⠀⠀⠙⢿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
+# ⠀⠀⠀⠀⠙⢿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⠀⠀⠀⠀⠙⢿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⠀⠀⠀⠀⠙⢿⣿⣿⣿⣿⣿⣿⣿⣿
+# ⠀⠀⠀⠀⠀⠂⠛⢿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⠀⠀⠀⠀⠀⠀⠙⢿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⠀⠀⠀⠀⠀⠀⠙⢿⣿⣿⣿⣿⣿⣿
+# ⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣿⣿⣿⣿⣿⣿⣿⣿⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣿⣿⣿⣿⣿⣿⣿⣿⠀⠀⠀⠀⠀⠀⢀⠐⠙⢿⣿⣿⣿⣿
+# ⠀⠀⠐⠀⠀⠀⠀⠀⠀⢀⢙⣿⣿⣿⣿⣿⣿⣿⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣿⣿⣿⣿⣿⣿⠀⠀⠀⠀⠀⠀⠁⠀⡀⠀⠙⢿⣿⣿
+# ⠀⠀⠐⠀⠀⠀⠀⠀⠀⠀⠄⠀⠙⢿⣿⣿⣿⣿⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣿⣿⣿⣿⠀⠀⠀⠀⠀⠀⡀⠀⠀⢀⠀⠀⠙⢿
+# ⠀⡀⠀⠀⠀⠀⠀⠀⠂⠒⠒⠀⠀⠀⠙⢿⣿⣿⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣿⣿⠀⠀⠀⠀⠀⠀⠀⢀⠀⠀⠀⡀⠀⠀
+# ⠀⠀⠀⠀⠀⠀⠀⠀⣈⡀⠀⠀⠀⠀⠀⠀⠙⢿⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠂⠀⢀⠀
+#
+function _triuppers_qr(T)
+    blocksize, cols = size(T)
+    return Iterators.map(0:div(cols - 1, blocksize)) do i
+        n = min(blocksize, cols - i * blocksize)
+        return UpperTriangular(view(T, 1:n, (1:n) .+ i * blocksize))
+    end
+end
+
 function Base.hash(F::QRCompactWY, h::UInt)
-    return hash(F.factors, hash(UpperTriangular(F.T), hash(QRCompactWY, h)))
+    return hash(F.factors, foldr(hash, _triuppers_qr(F.T); init=hash(QRCompactWY, h)))
 end
 function Base.:(==)(A::QRCompactWY, B::QRCompactWY)
-    return A.factors == B.factors && UpperTriangular(A.T) == UpperTriangular(B.T)
+    return A.factors == B.factors && all(splat(==), zip(_triuppers_qr.((A.T, B.T))...))
 end
 function Base.isequal(A::QRCompactWY, B::QRCompactWY)
-    return isequal(A.factors, B.factors) && isequal(UpperTriangular(A.T), UpperTriangular(B.T))
+    return isequal(A.factors, B.factors) && all(splat(isequal), zip(_triuppers_qr.((A.T, B.T))...))
 end
 
 """
diff --git a/stdlib/LinearAlgebra/test/factorization.jl b/stdlib/LinearAlgebra/test/factorization.jl
index a67871af4cdb0..6a9226d80cdf6 100644
--- a/stdlib/LinearAlgebra/test/factorization.jl
+++ b/stdlib/LinearAlgebra/test/factorization.jl
@@ -42,4 +42,19 @@ using Test, LinearAlgebra
     @test hash(F) == hash(G)
 end
 
+@testset "equality of QRCompactWY" begin
+    A = rand(100, 100)
+    F, G = qr(A), qr(A)
+
+    @test F == G
+    @test isequal(F, G)
+    @test hash(F) == hash(G)
+
+    G.T[28, 100] = 42
+
+    @test F != G
+    @test !isequal(F, G)
+    @test hash(F) != hash(G)
+end
+
 end

From 1c7dcbc8f4045e70cdfa8dc5e23d2a96a84c1097 Mon Sep 17 00:00:00 2001
From: Simeon Schaub <simeondavidschaub99@gmail.com>
Date: Mon, 28 Jun 2021 00:48:15 +0200
Subject: [PATCH 3/3] fix missings test

---
 stdlib/LinearAlgebra/src/qr.jl | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/stdlib/LinearAlgebra/src/qr.jl b/stdlib/LinearAlgebra/src/qr.jl
index 1ebf0231fb1ad..15bc61e1b1774 100644
--- a/stdlib/LinearAlgebra/src/qr.jl
+++ b/stdlib/LinearAlgebra/src/qr.jl
@@ -156,7 +156,9 @@ function Base.:(==)(A::QRCompactWY, B::QRCompactWY)
     return A.factors == B.factors && all(splat(==), zip(_triuppers_qr.((A.T, B.T))...))
 end
 function Base.isequal(A::QRCompactWY, B::QRCompactWY)
-    return isequal(A.factors, B.factors) && all(splat(isequal), zip(_triuppers_qr.((A.T, B.T))...))
+    return isequal(A.factors, B.factors) && all(zip(_triuppers_qr.((A.T, B.T))...)) do (a, b)
+        isequal(a, b)::Bool
+    end
 end
 
 """