From 4f99069df843d2f40918b05f0bfee482fc37198d Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreas@noack.dk>
Date: Wed, 19 May 2021 14:16:00 +0200
Subject: [PATCH 01/15] Add adjoint for Cholesky

---
 stdlib/LinearAlgebra/src/cholesky.jl  | 2 ++
 stdlib/LinearAlgebra/test/cholesky.jl | 8 ++++++++
 2 files changed, 10 insertions(+)

diff --git a/stdlib/LinearAlgebra/src/cholesky.jl b/stdlib/LinearAlgebra/src/cholesky.jl
index 18ee4cb5c7dd9..6e381243faf43 100644
--- a/stdlib/LinearAlgebra/src/cholesky.jl
+++ b/stdlib/LinearAlgebra/src/cholesky.jl
@@ -529,6 +529,8 @@ Base.propertynames(F::CholeskyPivoted, private::Bool=false) =
 
 issuccess(C::Union{Cholesky,CholeskyPivoted}) = C.info == 0
 
+adjoint(C::Union{Cholesky,CholeskyPivoted}) = C
+
 function show(io::IO, mime::MIME{Symbol("text/plain")}, C::Cholesky{<:Any,<:AbstractMatrix})
     if issuccess(C)
         summary(io, C); println(io)
diff --git a/stdlib/LinearAlgebra/test/cholesky.jl b/stdlib/LinearAlgebra/test/cholesky.jl
index a3f780c047a29..170af59eef8c3 100644
--- a/stdlib/LinearAlgebra/test/cholesky.jl
+++ b/stdlib/LinearAlgebra/test/cholesky.jl
@@ -475,4 +475,12 @@ end
     end
 end
 
+@testset "adjoint of Cholesky" begin
+    A = randn(5, 5)
+    A = A'A
+    F = cholesky(A)
+    b = ones(size(A, 1))
+    @test F\b == F'\b
+end
+
 end # module TestCholesky

From 077ef76abed4033dd1b8f08ff42ae81efb7d0b91 Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreas@noack.dk>
Date: Wed, 19 May 2021 14:29:47 +0200
Subject: [PATCH 02/15] Implement adjoint for BunchKaufman

---
 stdlib/LinearAlgebra/src/bunchkaufman.jl  | 10 +++++++++-
 stdlib/LinearAlgebra/test/bunchkaufman.jl | 21 ++++++++++++++++++++-
 2 files changed, 29 insertions(+), 2 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/bunchkaufman.jl b/stdlib/LinearAlgebra/src/bunchkaufman.jl
index 63254308799e1..daaa1f9dd61ea 100644
--- a/stdlib/LinearAlgebra/src/bunchkaufman.jl
+++ b/stdlib/LinearAlgebra/src/bunchkaufman.jl
@@ -205,7 +205,7 @@ convert(::Type{Factorization{T}}, B::BunchKaufman) where {T} = convert(BunchKauf
 size(B::BunchKaufman) = size(getfield(B, :LD))
 size(B::BunchKaufman, d::Integer) = size(getfield(B, :LD), d)
 issymmetric(B::BunchKaufman) = B.symmetric
-ishermitian(B::BunchKaufman) = !B.symmetric
+ishermitian(B::BunchKaufman{T}) where T = T<:Real || !B.symmetric
 
 function _ipiv2perm_bk(v::AbstractVector{T}, maxi::Integer, uplo::AbstractChar, rook::Bool) where T
     require_one_based_indexing(v)
@@ -279,6 +279,14 @@ Base.propertynames(B::BunchKaufman, private::Bool=false) =
 
 issuccess(B::BunchKaufman) = B.info == 0
 
+function adjoint(B::BunchKaufman)
+    if ishermitian(B)
+        return B
+    else
+        throw(ArgumentError("adjoint not implemented for complex symmetric matrices"))
+    end
+end
+
 function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, B::BunchKaufman)
     if issuccess(B)
         summary(io, B); println(io)
diff --git a/stdlib/LinearAlgebra/test/bunchkaufman.jl b/stdlib/LinearAlgebra/test/bunchkaufman.jl
index 5098f818f1804..e05592cdc1c5c 100644
--- a/stdlib/LinearAlgebra/test/bunchkaufman.jl
+++ b/stdlib/LinearAlgebra/test/bunchkaufman.jl
@@ -114,7 +114,8 @@ bimg  = randn(n,2)/2
                         @test logabsdet(bc2)[2] == sign(det(bc2))
                         @test inv(bc2)*apd ≈ Matrix(I, n, n)
                         @test apd*(bc2\b) ≈ b rtol=eps(cond(apd))
-                        @test ishermitian(bc2) == !issymmetric(bc2)
+                        @test ishermitian(bc2)
+                        @test !issymmetric(bc2) || eltya <: Real
                     end
                 end
             end
@@ -171,4 +172,22 @@ end
     end
 end
 
+@testset "adjoint of BunchKaufman" begin
+    Ar = randn(5, 5)
+    Ar = Ar + Ar'
+    Actmp = complex.(randn(5, 5), randn(5, 5))
+    Ac1 = Actmp + Actmp'
+    Ac2 = Actmp + transpose(Actmp)
+    b = ones(size(Ar, 1))
+
+    F = bunchkaufman(Ar)
+    @test F\b == F'\b
+
+    F = bunchkaufman(Ac1)
+    @test F\b == F'\b
+
+    F = bunchkaufman(Ac2)
+    @test_throws ArgumentError("adjoint not implemented for complex symmetric matrices") F'
+end
+
 end # module TestBunchKaufman

From 997cf8b259da3920c73e9f2d83380ea97e63cd6c Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreas@noack.dk>
Date: Wed, 19 May 2021 14:49:14 +0200
Subject: [PATCH 03/15] Fix ldiv! for adjoints of Hessenbergs

---
 stdlib/LinearAlgebra/src/hessenberg.jl  | 14 ++++++++------
 stdlib/LinearAlgebra/test/hessenberg.jl | 11 +++++++++++
 2 files changed, 19 insertions(+), 6 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/hessenberg.jl b/stdlib/LinearAlgebra/src/hessenberg.jl
index ffd12f39eafe2..e3b5d4983b30d 100644
--- a/stdlib/LinearAlgebra/src/hessenberg.jl
+++ b/stdlib/LinearAlgebra/src/hessenberg.jl
@@ -564,28 +564,30 @@ function AbstractMatrix(F::Hessenberg)
     end
 end
 
+# adjoint(Q::HessenbergQ{<:Real})
+
 lmul!(Q::BlasHessenbergQ{T,false}, X::StridedVecOrMat{T}) where {T<:BlasFloat} =
     LAPACK.ormhr!('L', 'N', 1, size(Q.factors, 1), Q.factors, Q.τ, X)
-rmul!(X::StridedMatrix{T}, Q::BlasHessenbergQ{T,false}) where {T<:BlasFloat} =
+rmul!(X::StridedVecOrMat{T}, Q::BlasHessenbergQ{T,false}) where {T<:BlasFloat} =
     LAPACK.ormhr!('R', 'N', 1, size(Q.factors, 1), Q.factors, Q.τ, X)
 lmul!(adjQ::Adjoint{<:Any,<:BlasHessenbergQ{T,false}}, X::StridedVecOrMat{T}) where {T<:BlasFloat} =
     (Q = adjQ.parent; LAPACK.ormhr!('L', ifelse(T<:Real, 'T', 'C'), 1, size(Q.factors, 1), Q.factors, Q.τ, X))
-rmul!(X::StridedMatrix{T}, adjQ::Adjoint{<:Any,<:BlasHessenbergQ{T,false}}) where {T<:BlasFloat} =
+rmul!(X::StridedVecOrMat{T}, adjQ::Adjoint{<:Any,<:BlasHessenbergQ{T,false}}) where {T<:BlasFloat} =
     (Q = adjQ.parent; LAPACK.ormhr!('R', ifelse(T<:Real, 'T', 'C'), 1, size(Q.factors, 1), Q.factors, Q.τ, X))
 
 lmul!(Q::BlasHessenbergQ{T,true}, X::StridedVecOrMat{T}) where {T<:BlasFloat} =
     LAPACK.ormtr!('L', Q.uplo, 'N', Q.factors, Q.τ, X)
-rmul!(X::StridedMatrix{T}, Q::BlasHessenbergQ{T,true}) where {T<:BlasFloat} =
+rmul!(X::StridedVecOrMat{T}, Q::BlasHessenbergQ{T,true}) where {T<:BlasFloat} =
     LAPACK.ormtr!('R', Q.uplo, 'N', Q.factors, Q.τ, X)
 lmul!(adjQ::Adjoint{<:Any,<:BlasHessenbergQ{T,true}}, X::StridedVecOrMat{T}) where {T<:BlasFloat} =
     (Q = adjQ.parent; LAPACK.ormtr!('L', Q.uplo, ifelse(T<:Real, 'T', 'C'), Q.factors, Q.τ, X))
-rmul!(X::StridedMatrix{T}, adjQ::Adjoint{<:Any,<:BlasHessenbergQ{T,true}}) where {T<:BlasFloat} =
+rmul!(X::StridedVecOrMat{T}, adjQ::Adjoint{<:Any,<:BlasHessenbergQ{T,true}}) where {T<:BlasFloat} =
     (Q = adjQ.parent; LAPACK.ormtr!('R', Q.uplo, ifelse(T<:Real, 'T', 'C'), Q.factors, Q.τ, X))
 
 lmul!(Q::HessenbergQ{T}, X::Adjoint{T,<:StridedVecOrMat{T}}) where {T} = rmul!(X', Q')'
-rmul!(X::Adjoint{T,<:StridedMatrix{T}}, Q::HessenbergQ{T}) where {T} = lmul!(Q', X')'
+rmul!(X::Adjoint{T,<:StridedVecOrMat{T}}, Q::HessenbergQ{T}) where {T} = lmul!(Q', X')'
 lmul!(adjQ::Adjoint{<:Any,<:HessenbergQ{T}}, X::Adjoint{T,<:StridedVecOrMat{T}}) where {T}  = rmul!(X', adjQ')'
-rmul!(X::Adjoint{T,<:StridedMatrix{T}}, adjQ::Adjoint{<:Any,<:HessenbergQ{T}}) where {T} = lmul!(adjQ', X')'
+rmul!(X::Adjoint{T,<:StridedVecOrMat{T}}, adjQ::Adjoint{<:Any,<:HessenbergQ{T}}) where {T} = lmul!(adjQ', X')'
 
 # multiply x by the entries of M in the upper-k triangle, which contains
 # the entries of the upper-Hessenberg matrix H for k=-1
diff --git a/stdlib/LinearAlgebra/test/hessenberg.jl b/stdlib/LinearAlgebra/test/hessenberg.jl
index dd6a131a6ae1e..0f9246c722349 100644
--- a/stdlib/LinearAlgebra/test/hessenberg.jl
+++ b/stdlib/LinearAlgebra/test/hessenberg.jl
@@ -185,4 +185,15 @@ end
     @test Base.propertynames(F, true) == (:Q, :H, :μ, :τ, :factors, :uplo)
 end
 
+@testset "adjoint of Hessenberg" begin
+    Ar = randn(5, 5)
+    Ac = complex.(randn(5, 5), randn(5, 5))
+    b = ones(size(Ar, 1))
+
+    for A in (Ar, Ac)
+        F = hessenberg(A)
+        @test A'\b ≈ F'\b
+    end
+end
+
 end # module TestHessenberg

From f815e99e8483f73b6f4e36281eaf562a3129b540 Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreas@noack.dk>
Date: Wed, 19 May 2021 15:04:08 +0200
Subject: [PATCH 04/15] Add adjoint of LDLt

---
 stdlib/LinearAlgebra/src/ldlt.jl     |  3 +++
 stdlib/LinearAlgebra/test/tridiag.jl | 12 ++++++++++++
 2 files changed, 15 insertions(+)

diff --git a/stdlib/LinearAlgebra/src/ldlt.jl b/stdlib/LinearAlgebra/src/ldlt.jl
index d0f59ebb9ff1b..f1ea10aa0f614 100644
--- a/stdlib/LinearAlgebra/src/ldlt.jl
+++ b/stdlib/LinearAlgebra/src/ldlt.jl
@@ -77,6 +77,9 @@ function getproperty(F::LDLt, d::Symbol)
     end
 end
 
+adjoint(F::LDLt{<:Real,<:SymTridiagonal}) = F
+adjoint(F::LDLt) = LDLt(copy(adjoint(F.data)))
+
 function show(io::IO, mime::MIME{Symbol("text/plain")}, F::LDLt)
     summary(io, F); println(io)
     println(io, "L factor:")
diff --git a/stdlib/LinearAlgebra/test/tridiag.jl b/stdlib/LinearAlgebra/test/tridiag.jl
index 8cfde8fb5067c..d4137e3738206 100644
--- a/stdlib/LinearAlgebra/test/tridiag.jl
+++ b/stdlib/LinearAlgebra/test/tridiag.jl
@@ -605,4 +605,16 @@ end
     end
 end
 
+@testset "adjoint of LDLt" begin
+    Sr = SymTridiagonal(randn(5), randn(4))
+    Sc = SymTridiagonal(complex.(randn(5)) .+ 1im, complex.(randn(4), randn(4)))
+    b = ones(size(Sr, 1))
+
+    F = ldlt(Sr)
+    @test F\b == F'\b
+
+    F = ldlt(Sc)
+    @test copy(Sc')\b == F'\b
+end
+
 end # module TestTridiagonal

From a3755d8090615d99564a79c0e2b077694e86cf61 Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreas@noack.dk>
Date: Thu, 20 May 2021 21:09:32 +0200
Subject: [PATCH 05/15] Fix return for tall problems in fallback \ method for
 adjoint of Factorizations to make \ work for adjoint LQ.

---
 stdlib/LinearAlgebra/src/factorization.jl |  6 ++++++
 stdlib/LinearAlgebra/src/lq.jl            | 26 ++++++++++++-----------
 stdlib/LinearAlgebra/test/lq.jl           | 19 +++++++++++++++++
 3 files changed, 39 insertions(+), 12 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/factorization.jl b/stdlib/LinearAlgebra/src/factorization.jl
index 5ff215a3eb665..79e67d4dc0476 100644
--- a/stdlib/LinearAlgebra/src/factorization.jl
+++ b/stdlib/LinearAlgebra/src/factorization.jl
@@ -103,6 +103,8 @@ function \(F::Factorization, B::AbstractVecOrMat)
     copyto!(BB, B)
     ldiv!(F, BB)
 end
+_rows(b::AbstractVector, r::AbstractVector) = b[r]
+_rows(B::AbstractVector, r::AbstractMatrix) = B[r, :]
 function \(adjF::Adjoint{<:Any,<:Factorization}, B::AbstractVecOrMat)
     require_one_based_indexing(B)
     F = adjF.parent
@@ -110,6 +112,10 @@ function \(adjF::Adjoint{<:Any,<:Factorization}, B::AbstractVecOrMat)
     BB = similar(B, TFB, size(B))
     copyto!(BB, B)
     ldiv!(adjoint(F), BB)
+    # For tall problems, we compute a least sqaures solution so only part
+    # of the rhs should be returned from \ while ldiv! uses (and returns)
+    # the complete rhs
+    return (>)(size(adjF)...) ? _rows(BB, 1:size(adjF, 2)) : BB
 end
 
 function /(B::AbstractMatrix, F::Factorization)
diff --git a/stdlib/LinearAlgebra/src/lq.jl b/stdlib/LinearAlgebra/src/lq.jl
index 606e8c3dd006e..f586e374a542d 100644
--- a/stdlib/LinearAlgebra/src/lq.jl
+++ b/stdlib/LinearAlgebra/src/lq.jl
@@ -318,17 +318,6 @@ _rightappdimmismatch(rowsorcols) =
         "or (2) the number of rows of that (LQPackedQ) matrix's internal representation ",
         "(the factorization's originating matrix's number of rows)")))
 
-
-function (\)(A::LQ{TA},B::StridedVecOrMat{TB}) where {TA,TB}
-    S = promote_type(TA,TB)
-    m, n = size(A)
-    m ≤ n || throw(DimensionMismatch("LQ solver does not support overdetermined systems (more rows than columns)"))
-    m == size(B,1) || throw(DimensionMismatch("Both inputs should have the same number of rows"))
-    AA = Factorization{S}(A)
-    X = _zeros(S, B, n)
-    X[1:size(B, 1), :] = B
-    return ldiv!(AA, X)
-end
 # With a real lhs and complex rhs with the same precision, we can reinterpret
 # the complex rhs as a real rhs with twice the number of columns
 function (\)(F::LQ{T}, B::VecOrMat{Complex{T}}) where T<:BlasReal
@@ -342,12 +331,25 @@ function (\)(F::LQ{T}, B::VecOrMat{Complex{T}}) where T<:BlasReal
 end
 
 
-function ldiv!(A::LQ{T}, B::StridedVecOrMat{T}) where T
+function ldiv!(A::LQ, B::StridedVecOrMat)
     require_one_based_indexing(B)
+    m, n = size(A)
+    m ≤ n || throw(DimensionMismatch("LQ solver does not support overdetermined systems (more rows than columns)"))
+
     ldiv!(LowerTriangular(A.L), view(B, 1:size(A,1), axes(B,2)))
     return lmul!(adjoint(A.Q), B)
 end
 
+function ldiv!(Fadj::Adjoint{<:Any,<:LQ}, B::StridedVecOrMat)
+    require_one_based_indexing(B)
+    m, n = size(Fadj)
+    m >= n || throw(DimensionMismatch("solver does not support underdetermined systems (more columns than rows)"))
+
+    F = parent(Fadj)
+    lmul!(F.Q, B)
+    ldiv!(UpperTriangular(adjoint(F.L)), view(B, 1:size(F,1), axes(B,2)))
+    return B
+end
 
 # In LQ factorization, `Q` is expressed as the product of the adjoint of the
 # reflectors.  Thus, `det` has to be conjugated.
diff --git a/stdlib/LinearAlgebra/test/lq.jl b/stdlib/LinearAlgebra/test/lq.jl
index 8915324af9461..1cf4a26f906f6 100644
--- a/stdlib/LinearAlgebra/test/lq.jl
+++ b/stdlib/LinearAlgebra/test/lq.jl
@@ -220,4 +220,23 @@ Q factor:
  0.0  0.0  0.0  1.0"""
 end
 
+@testset "adjoint of LQ" begin
+    n = 5
+    B = ones(n, 2)
+
+    for b in (B[:, 1], B)
+        for A in (
+            randn(n, n),
+            # Tall problems become least squares problems similarly to QR
+            randn(n - 2, n),
+            complex.(randn(n, n), randn(n, n)))
+
+            F = lq(A)
+            @test A'\b ≈ F'\b
+        end
+        @test_throws DimensionMismatch lq(randn(n, n + 2))'\b
+    end
+
+end
+
 end # module TestLQ

From c0446089bd7ead5d647ffe8b54c606db714c8677 Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreas@noack.dk>
Date: Fri, 21 May 2021 10:50:49 +0200
Subject: [PATCH 06/15] Fix qr(A)'\b

---
 stdlib/LinearAlgebra/src/factorization.jl | 50 +++++++++++++++--------
 stdlib/LinearAlgebra/src/qr.jl            | 45 ++++++++++++--------
 stdlib/LinearAlgebra/test/qr.jl           | 44 ++++++++++++++++++--
 3 files changed, 102 insertions(+), 37 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/factorization.jl b/stdlib/LinearAlgebra/src/factorization.jl
index 79e67d4dc0476..d4009e29047bb 100644
--- a/stdlib/LinearAlgebra/src/factorization.jl
+++ b/stdlib/LinearAlgebra/src/factorization.jl
@@ -96,26 +96,44 @@ function (/)(B::VecOrMat{Complex{T}}, F::Factorization{T}) where T<:BlasReal
     return copy(reinterpret(Complex{T}, x))
 end
 
-function \(F::Factorization, B::AbstractVecOrMat)
+# convenience methods
+## return only the solution of a least squares problem while avoiding promoting
+## vectors to matrices.
+_cut_B(x::AbstractVector, r::UnitRange) = length(x)  > length(r) ? x[r]   : x
+_cut_B(X::AbstractMatrix, r::UnitRange) = size(X, 1) > length(r) ? X[r,:] : X
+
+## append right hand side with zeros if necessary
+_zeros(::Type{T}, b::AbstractVector, n::Integer) where {T} = zeros(T, max(length(b), n))
+_zeros(::Type{T}, B::AbstractMatrix, n::Integer) where {T} = zeros(T, max(size(B, 1), n), size(B, 2))
+
+# General fallback definition for handling under- and overdetermined system as well as square problems
+function \(adjF::Union{<:Factorization,Adjoint{<:Any,<:Factorization}}, B::AbstractVecOrMat)
     require_one_based_indexing(B)
-    TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(F)))
-    BB = similar(B, TFB, size(B))
-    copyto!(BB, B)
-    ldiv!(F, BB)
-end
-_rows(b::AbstractVector, r::AbstractVector) = b[r]
-_rows(B::AbstractVector, r::AbstractMatrix) = B[r, :]
-function \(adjF::Adjoint{<:Any,<:Factorization}, B::AbstractVecOrMat)
-    require_one_based_indexing(B)
-    F = adjF.parent
-    TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(F)))
-    BB = similar(B, TFB, size(B))
-    copyto!(BB, B)
-    ldiv!(adjoint(F), BB)
+    m, n = size(adjF)
+    if m != size(B, 1)
+        throw(DimensionMismatch("arguments must have the same number of rows"))
+    end
+    # F = adjF.parent
+    TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(adjF)))
+
+    # For wide problem we (often) compute a minimum norm solution. The solution
+    # is larger than the right hand side so we use size(adjF, 2).
+    BB = _zeros(TFB, B, n)
+
+    if n > size(B, 1)
+        # Underdetermined
+        fill!(BB, 0)
+        copyto!(view(BB, 1:m, :), B)
+    else
+        copyto!(BB, B)
+    end
+
+    ldiv!(adjF, BB)
+
     # For tall problems, we compute a least sqaures solution so only part
     # of the rhs should be returned from \ while ldiv! uses (and returns)
     # the complete rhs
-    return (>)(size(adjF)...) ? _rows(BB, 1:size(adjF, 2)) : BB
+    return _cut_B(BB, 1:n)
 end
 
 function /(B::AbstractMatrix, F::Factorization)
diff --git a/stdlib/LinearAlgebra/src/qr.jl b/stdlib/LinearAlgebra/src/qr.jl
index 390c8a5875773..a2f8a49639c9e 100644
--- a/stdlib/LinearAlgebra/src/qr.jl
+++ b/stdlib/LinearAlgebra/src/qr.jl
@@ -472,6 +472,8 @@ end
 Base.propertynames(F::QRPivoted, private::Bool=false) =
     (:R, :Q, :p, :P, (private ? fieldnames(typeof(F)) : ())...)
 
+adjoint(F::Union{QR,QRPivoted,QRCompactWY}) = Adjoint(F)
+
 abstract type AbstractQ{T} <: AbstractMatrix{T} end
 
 inv(Q::AbstractQ) = Q'
@@ -939,28 +941,35 @@ function ldiv!(A::QRPivoted, B::StridedMatrix)
     B
 end
 
-# convenience methods
-## return only the solution of a least squares problem while avoiding promoting
-## vectors to matrices.
-_cut_B(x::AbstractVector, r::UnitRange) = length(x)  > length(r) ? x[r]   : x
-_cut_B(X::AbstractMatrix, r::UnitRange) = size(X, 1) > length(r) ? X[r,:] : X
-
-## append right hand side with zeros if necessary
-_zeros(::Type{T}, b::AbstractVector, n::Integer) where {T} = zeros(T, max(length(b), n))
-_zeros(::Type{T}, B::AbstractMatrix, n::Integer) where {T} = zeros(T, max(size(B, 1), n), size(B, 2))
+function _apply_permutation!(F::QRPivoted, B::AbstractVecOrMat)
+    # Apply permutation but only to the top part of the solution vector since
+    # it's padded with zeros for underdetermined problems
+    B[1:length(F.p), :] = B[F.p, :]
+    return B
+end
+_apply_permutation!(F::Factorization, B::AbstractVecOrMat) = B
 
-function (\)(A::Union{QR{TA},QRCompactWY{TA},QRPivoted{TA}}, B::AbstractVecOrMat{TB}) where {TA,TB}
+function ldiv!(Fadj::Adjoint{<:Any,<:Union{QR,QRCompactWY,QRPivoted}}, B::AbstractVecOrMat)
     require_one_based_indexing(B)
-    S = promote_type(TA,TB)
-    m, n = size(A)
-    m == size(B,1) || throw(DimensionMismatch("Both inputs should have the same number of rows"))
+    m, n = size(Fadj)
 
-    AA = Factorization{S}(A)
+    # We don't allow solutions overdetermined systems. It would at least be
+    if m > n
+        throw(DimensionMismatch("overdetermined systems are not supported"))
+    end
+    if n != size(B, 1)
+        throw(DimensionMismatch("inputs should have the same number of rows"))
+    end
+    F = parent(Fadj)
 
-    X = _zeros(S, B, n)
-    X[1:size(B, 1), :] = B
-    ldiv!(AA, X)
-    return _cut_B(X, 1:n)
+    B = _apply_permutation!(F, B)
+
+    # For underdetermined system, the triangular solve should only be applied to the top
+    # part of B that contains the rhs. For square problems, the view corresponds to B itself
+    ldiv!(LowerTriangular(adjoint(F.R)), view(B, 1:size(F.R, 2), :))
+    lmul!(F.Q, B)
+
+    return B
 end
 
 # With a real lhs and complex rhs with the same precision, we can reinterpret the complex
diff --git a/stdlib/LinearAlgebra/test/qr.jl b/stdlib/LinearAlgebra/test/qr.jl
index 16f828b4f8861..b713d01ab4e8d 100644
--- a/stdlib/LinearAlgebra/test/qr.jl
+++ b/stdlib/LinearAlgebra/test/qr.jl
@@ -212,11 +212,8 @@ end
 
 @testset "transpose errors" begin
     @test_throws MethodError transpose(qr(randn(3,3)))
-    @test_throws MethodError adjoint(qr(randn(3,3)))
     @test_throws MethodError transpose(qr(randn(3,3), NoPivot()))
-    @test_throws MethodError adjoint(qr(randn(3,3), NoPivot()))
     @test_throws MethodError transpose(qr(big.(randn(3,3))))
-    @test_throws MethodError adjoint(qr(big.(randn(3,3))))
 end
 
 @testset "Issue 7304" begin
@@ -369,4 +366,45 @@ end
     end
 end
 
+@testset "adjoint of QR" begin
+    n = 5
+    B = randn(5, 2)
+
+    @testset "size(b)=$(size(b))" for b in (B[:, 1], B)
+        @testset "size(A)=$(size(A))" for A in (
+            randn(n, n),
+            # Wide problems become minimum norm (in x) problems similarly to LQ
+            randn(n + 2, n),
+            complex.(randn(n, n), randn(n, n)))
+
+            @testset "QRCompactWY" begin
+                F = qr(A)
+                x = F'\b
+                @test x ≈ A'\b
+                @test length(size(x)) == length(size(b))
+            end
+
+            @testset "QR" begin
+                F = LinearAlgebra.qrfactUnblocked!(copy(A))
+                x = F'\b
+                @test x ≈ A'\b
+                @test length(size(x)) == length(size(b))
+            end
+
+            @testset "QRPivoted" begin
+                F = LinearAlgebra.qr(A, Val(true))
+                x = F'\b
+                @test x ≈ A'\b
+                @test length(size(x)) == length(size(b))
+            end
+        end
+        @test_throws DimensionMismatch("overdetermined systems are not supported")    qr(randn(n - 2, n))'\b
+        @test_throws DimensionMismatch("arguments must have the same number of rows") qr(randn(n, n + 1))'\b
+        @test_throws DimensionMismatch("overdetermined systems are not supported")    LinearAlgebra.qrfactUnblocked!(randn(n - 2, n))'\b
+        @test_throws DimensionMismatch("arguments must have the same number of rows") LinearAlgebra.qrfactUnblocked!(randn(n, n + 1))'\b
+        @test_throws DimensionMismatch("overdetermined systems are not supported")    qr(randn(n - 2, n), Val(true))'\b
+        @test_throws DimensionMismatch("arguments must have the same number of rows") qr(randn(n, n + 1), Val(true))'\b
+    end
+end
+
 end # module TestQR

From 8dd97817a9cc1d79239a3613a2b486feaf73f870 Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreas@noack.dk>
Date: Fri, 21 May 2021 14:30:50 +0200
Subject: [PATCH 07/15] Define adjoint for SVD

---
 stdlib/LinearAlgebra/src/svd.jl  |  4 ++++
 stdlib/LinearAlgebra/test/svd.jl | 20 ++++++++++++++++++++
 2 files changed, 24 insertions(+)

diff --git a/stdlib/LinearAlgebra/src/svd.jl b/stdlib/LinearAlgebra/src/svd.jl
index 68bce4793661f..c5063a6117547 100644
--- a/stdlib/LinearAlgebra/src/svd.jl
+++ b/stdlib/LinearAlgebra/src/svd.jl
@@ -252,6 +252,10 @@ end
 size(A::SVD, dim::Integer) = dim == 1 ? size(A.U, dim) : size(A.Vt, dim)
 size(A::SVD) = (size(A, 1), size(A, 2))
 
+function adjoint(F::SVD)
+    return SVD(F.Vt', F.S, F.U')
+end
+
 function show(io::IO, mime::MIME{Symbol("text/plain")}, F::SVD{<:Any,<:Any,<:AbstractArray})
     summary(io, F); println(io)
     println(io, "U factor:")
diff --git a/stdlib/LinearAlgebra/test/svd.jl b/stdlib/LinearAlgebra/test/svd.jl
index 30dd6db300eb9..f02b8def49e82 100644
--- a/stdlib/LinearAlgebra/test/svd.jl
+++ b/stdlib/LinearAlgebra/test/svd.jl
@@ -217,4 +217,24 @@ end
     @test Uc * diagm(0=>Sc) * transpose(V) ≈ complex.(A) rtol=1e-3
 end
 
+@testset "adjoint of SVD" begin
+    n = 5
+    B = randn(5, 2)
+
+    @testset "size(b)=$(size(b))" for b in (B[:, 1], B)
+        @testset "size(A)=$(size(A))" for A in (
+            randn(n, n),
+            # Wide problems become minimum norm (in x) problems similarly to LQ
+            randn(n + 2, n),
+            randn(n - 2, n),
+            complex.(randn(n, n), randn(n, n)))
+
+            F = svd(A)
+            x = F'\b
+            @test x ≈ A'\b
+            @test length(size(x)) == length(size(b))
+        end
+    end
+end
+
 end # module TestSVD

From a350a7627ad5fc5e22d04b22991d6eaf4ece1233 Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreas@noack.dk>
Date: Mon, 24 May 2021 15:32:05 +0200
Subject: [PATCH 08/15] Improve promotion in fallback by defining general
 convert methods for Factorizations

---
 stdlib/LinearAlgebra/src/bunchkaufman.jl  |  6 ++----
 stdlib/LinearAlgebra/src/factorization.jl | 16 ++++++++++------
 stdlib/LinearAlgebra/src/lq.jl            |  4 ++--
 stdlib/LinearAlgebra/test/lq.jl           |  3 +--
 4 files changed, 15 insertions(+), 14 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/bunchkaufman.jl b/stdlib/LinearAlgebra/src/bunchkaufman.jl
index daaa1f9dd61ea..75fb9ae7bf04e 100644
--- a/stdlib/LinearAlgebra/src/bunchkaufman.jl
+++ b/stdlib/LinearAlgebra/src/bunchkaufman.jl
@@ -196,11 +196,9 @@ julia> S.L*S.D*S.L' - A[S.p, S.p]
 bunchkaufman(A::AbstractMatrix{T}, rook::Bool=false; check::Bool = true) where {T} =
     bunchkaufman!(copy_oftype(A, typeof(sqrt(oneunit(T)))), rook; check = check)
 
-convert(::Type{BunchKaufman{T}}, B::BunchKaufman{T}) where {T} = B
-convert(::Type{BunchKaufman{T}}, B::BunchKaufman) where {T} =
+BunchKaufman{T}(B::BunchKaufman) where {T} =
     BunchKaufman(convert(Matrix{T}, B.LD), B.ipiv, B.uplo, B.symmetric, B.rook, B.info)
-convert(::Type{Factorization{T}}, B::BunchKaufman{T}) where {T} = B
-convert(::Type{Factorization{T}}, B::BunchKaufman) where {T} = convert(BunchKaufman{T}, B)
+Factorization{T}(B::BunchKaufman) where {T} = BunchKaufman{T}(B)
 
 size(B::BunchKaufman) = size(getfield(B, :LD))
 size(B::BunchKaufman, d::Integer) = size(getfield(B, :LD), d)
diff --git a/stdlib/LinearAlgebra/src/factorization.jl b/stdlib/LinearAlgebra/src/factorization.jl
index d4009e29047bb..a5d7873d27eea 100644
--- a/stdlib/LinearAlgebra/src/factorization.jl
+++ b/stdlib/LinearAlgebra/src/factorization.jl
@@ -59,6 +59,9 @@ convert(::Type{T}, f::Factorization) where {T<:AbstractArray} = T(f)
 
 ### General promotion rules
 Factorization{T}(F::Factorization{T}) where {T} = F
+# This is a bit odd since the return is not a Factorization but it works well in generic code
+Factorization{T}(A::Adjoint{<:Any,<:Factorization}) where {T} =
+    adjoint(Factorization{T}(parent(A)))
 inv(F::Factorization{T}) where {T} = (n = size(F, 1); ldiv!(F, Matrix{T}(I, n, n)))
 
 Base.hash(F::Factorization, h::UInt) = mapreduce(f -> hash(getfield(F, f)), hash, 1:nfields(F); init=h)
@@ -107,17 +110,18 @@ _zeros(::Type{T}, b::AbstractVector, n::Integer) where {T} = zeros(T, max(length
 _zeros(::Type{T}, B::AbstractMatrix, n::Integer) where {T} = zeros(T, max(size(B, 1), n), size(B, 2))
 
 # General fallback definition for handling under- and overdetermined system as well as square problems
-function \(adjF::Union{<:Factorization,Adjoint{<:Any,<:Factorization}}, B::AbstractVecOrMat)
+function \(F::Union{<:Factorization,Adjoint{<:Any,<:Factorization}}, B::AbstractVecOrMat)
     require_one_based_indexing(B)
-    m, n = size(adjF)
+    m, n = size(F)
     if m != size(B, 1)
         throw(DimensionMismatch("arguments must have the same number of rows"))
     end
-    # F = adjF.parent
-    TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(adjF)))
+
+    TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(F)))
+    FF = Factorization{TFB}(F)
 
     # For wide problem we (often) compute a minimum norm solution. The solution
-    # is larger than the right hand side so we use size(adjF, 2).
+    # is larger than the right hand side so we use size(F, 2).
     BB = _zeros(TFB, B, n)
 
     if n > size(B, 1)
@@ -128,7 +132,7 @@ function \(adjF::Union{<:Factorization,Adjoint{<:Any,<:Factorization}}, B::Abstr
         copyto!(BB, B)
     end
 
-    ldiv!(adjF, BB)
+    ldiv!(FF, BB)
 
     # For tall problems, we compute a least sqaures solution so only part
     # of the rhs should be returned from \ while ldiv! uses (and returns)
diff --git a/stdlib/LinearAlgebra/src/lq.jl b/stdlib/LinearAlgebra/src/lq.jl
index f586e374a542d..d9e49f3a9c1c6 100644
--- a/stdlib/LinearAlgebra/src/lq.jl
+++ b/stdlib/LinearAlgebra/src/lq.jl
@@ -125,8 +125,8 @@ lq_eltype(::Type{T}) where {T} = typeof(zero(T) / sqrt(abs2(one(T))))
 copy(A::LQ) = LQ(copy(A.factors), copy(A.τ))
 
 LQ{T}(A::LQ) where {T} = LQ(convert(AbstractMatrix{T}, A.factors), convert(Vector{T}, A.τ))
-Factorization{T}(A::LQ{T}) where {T} = A
 Factorization{T}(A::LQ) where {T} = LQ{T}(A)
+
 AbstractMatrix(A::LQ) = A.L*A.Q
 AbstractArray(A::LQ) = AbstractMatrix(A)
 Matrix(A::LQ) = Array(AbstractArray(A))
@@ -194,7 +194,7 @@ function lmul!(A::LQ, B::StridedVecOrMat)
 end
 function *(A::LQ{TA}, B::StridedVecOrMat{TB}) where {TA,TB}
     TAB = promote_type(TA, TB)
-    _cut_B(lmul!(Factorization{TAB}(A), copy_oftype(B, TAB)), 1:size(A,1))
+    _cut_B(lmul!(convert(Factorization{TAB}, A), copy_oftype(B, TAB)), 1:size(A,1))
 end
 
 ## Multiplication by Q
diff --git a/stdlib/LinearAlgebra/test/lq.jl b/stdlib/LinearAlgebra/test/lq.jl
index 1cf4a26f906f6..b054621e11313 100644
--- a/stdlib/LinearAlgebra/test/lq.jl
+++ b/stdlib/LinearAlgebra/test/lq.jl
@@ -222,9 +222,8 @@ end
 
 @testset "adjoint of LQ" begin
     n = 5
-    B = ones(n, 2)
 
-    for b in (B[:, 1], B)
+    for b in (ones(n), ones(n, 2), ones(Complex{Float64}, n, 2))
         for A in (
             randn(n, n),
             # Tall problems become least squares problems similarly to QR

From a1095581af679eedcd30431872033fc6e137aef8 Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreas@noack.dk>
Date: Mon, 24 May 2021 22:06:54 +0200
Subject: [PATCH 09/15] Fix ldiv! for SVD

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

diff --git a/stdlib/LinearAlgebra/src/svd.jl b/stdlib/LinearAlgebra/src/svd.jl
index c5063a6117547..6991cc9bf13d4 100644
--- a/stdlib/LinearAlgebra/src/svd.jl
+++ b/stdlib/LinearAlgebra/src/svd.jl
@@ -237,8 +237,9 @@ svdvals(S::SVD{<:Any,T}) where {T} = (S.S)::Vector{T}
 
 # SVD least squares
 function ldiv!(A::SVD{T}, B::StridedVecOrMat) where T
+    m, n = size(A)
     k = searchsortedlast(A.S, eps(real(T))*A.S[1], rev=true)
-    view(A.Vt,1:k,:)' * (view(A.S,1:k) .\ (view(A.U,:,1:k)' * B))
+    return mul!(view(B, 1:n, :), view(A.Vt, 1:k, :)', view(A.S, 1:k) .\ (view(A.U, :, 1:k)' * _cut_B(B, 1:m)))
 end
 
 function inv(F::SVD{T}) where T

From 0245e6c2d2c1d171eff1d02c91b769630f6e2cff Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreas@noack.dk>
Date: Tue, 25 May 2021 16:54:28 +0200
Subject: [PATCH 10/15] Restrict the general \ definition that handles over-
 and underdetermined systems to LAPACK factorizations

---
 stdlib/LinearAlgebra/src/LinearAlgebra.jl | 55 +++++++++++++++++++++++
 stdlib/LinearAlgebra/src/factorization.jl | 52 +++------------------
 stdlib/LinearAlgebra/src/svd.jl           |  2 +-
 3 files changed, 62 insertions(+), 47 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl
index d5a2a64467f93..98638fc6de7cf 100644
--- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl
+++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl
@@ -397,6 +397,61 @@ const ⋅ = dot
 const × = cross
 export ⋅, ×
 
+## convenience methods
+## return only the solution of a least squares problem while avoiding promoting
+## vectors to matrices.
+_cut_B(x::AbstractVector, r::UnitRange) = length(x)  > length(r) ? x[r]   : x
+_cut_B(X::AbstractMatrix, r::UnitRange) = size(X, 1) > length(r) ? X[r,:] : X
+
+## append right hand side with zeros if necessary
+_zeros(::Type{T}, b::AbstractVector, n::Integer) where {T} = zeros(T, max(length(b), n))
+_zeros(::Type{T}, B::AbstractMatrix, n::Integer) where {T} = zeros(T, max(size(B, 1), n), size(B, 2))
+
+# General fallback definition for handling under- and overdetermined system as well as square problems
+# While this definition is pretty general, it does e.g. promote to common element type of lhs and rhs
+# which is required by LAPACK but not SuiteSpase which allows real-complex solves in some cases. Hence,
+# we restrict this method to only the LAPACK factorizations in LinearAlgebra.
+# The definition is put here since it explicitly references all the Factorizion structs so it has
+# to be located after all the files that define the structs.
+const LAPACKFactorizations{T,S} = Union{
+    BunchKaufman{T,S},
+    Cholesky{T,S},
+    LQ{T,S},
+    LU{T,S},
+    QR{T,S},
+    QRCompactWY{T,S},
+    QRPivoted{T,S},
+    SVD{T,<:Real,S}}
+function (\)(F::Union{<:LAPACKFactorizations,Adjoint{<:Any,<:LAPACKFactorizations}}, B::AbstractVecOrMat)
+    require_one_based_indexing(B)
+    m, n = size(F)
+    if m != size(B, 1)
+        throw(DimensionMismatch("arguments must have the same number of rows"))
+    end
+
+    TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(F)))
+    FF = Factorization{TFB}(F)
+
+    # For wide problem we (often) compute a minimum norm solution. The solution
+    # is larger than the right hand side so we use size(F, 2).
+    BB = _zeros(TFB, B, n)
+
+    if n > size(B, 1)
+        # Underdetermined
+        fill!(BB, 0)
+        copyto!(view(BB, 1:m, :), B)
+    else
+        copyto!(BB, B)
+    end
+
+    ldiv!(FF, BB)
+
+    # For tall problems, we compute a least sqaures solution so only part
+    # of the rhs should be returned from \ while ldiv! uses (and returns)
+    # the complete rhs
+    return _cut_B(BB, 1:n)
+end
+
 """
     LinearAlgebra.peakflops(n::Integer=2000; parallel::Bool=false)
 
diff --git a/stdlib/LinearAlgebra/src/factorization.jl b/stdlib/LinearAlgebra/src/factorization.jl
index a5d7873d27eea..b651e85512f6d 100644
--- a/stdlib/LinearAlgebra/src/factorization.jl
+++ b/stdlib/LinearAlgebra/src/factorization.jl
@@ -99,65 +99,25 @@ function (/)(B::VecOrMat{Complex{T}}, F::Factorization{T}) where T<:BlasReal
     return copy(reinterpret(Complex{T}, x))
 end
 
-# convenience methods
-## return only the solution of a least squares problem while avoiding promoting
-## vectors to matrices.
-_cut_B(x::AbstractVector, r::UnitRange) = length(x)  > length(r) ? x[r]   : x
-_cut_B(X::AbstractMatrix, r::UnitRange) = size(X, 1) > length(r) ? X[r,:] : X
-
-## append right hand side with zeros if necessary
-_zeros(::Type{T}, b::AbstractVector, n::Integer) where {T} = zeros(T, max(length(b), n))
-_zeros(::Type{T}, B::AbstractMatrix, n::Integer) where {T} = zeros(T, max(size(B, 1), n), size(B, 2))
-
-# General fallback definition for handling under- and overdetermined system as well as square problems
-function \(F::Union{<:Factorization,Adjoint{<:Any,<:Factorization}}, B::AbstractVecOrMat)
-    require_one_based_indexing(B)
-    m, n = size(F)
-    if m != size(B, 1)
-        throw(DimensionMismatch("arguments must have the same number of rows"))
-    end
-
-    TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(F)))
-    FF = Factorization{TFB}(F)
-
-    # For wide problem we (often) compute a minimum norm solution. The solution
-    # is larger than the right hand side so we use size(F, 2).
-    BB = _zeros(TFB, B, n)
-
-    if n > size(B, 1)
-        # Underdetermined
-        fill!(BB, 0)
-        copyto!(view(BB, 1:m, :), B)
-    else
-        copyto!(BB, B)
-    end
-
-    ldiv!(FF, BB)
-
-    # For tall problems, we compute a least sqaures solution so only part
-    # of the rhs should be returned from \ while ldiv! uses (and returns)
-    # the complete rhs
-    return _cut_B(BB, 1:n)
-end
-
-function /(B::AbstractMatrix, F::Factorization)
+function \(F::Union{Factorization, Adjoint{<:Any,<:Factorization}}, B::AbstractVecOrMat)
     require_one_based_indexing(B)
     TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(F)))
     BB = similar(B, TFB, size(B))
     copyto!(BB, B)
-    rdiv!(BB, F)
+    ldiv!(F, BB)
 end
-function /(B::AbstractMatrix, adjF::Adjoint{<:Any,<:Factorization})
+
+function /(B::AbstractMatrix, F::Union{Factorization, Adjoint{<:Any,<:Factorization}})
     require_one_based_indexing(B)
-    F = adjF.parent
     TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(F)))
     BB = similar(B, TFB, size(B))
     copyto!(BB, B)
-    rdiv!(BB, adjoint(F))
+    rdiv!(BB, F)
 end
 /(adjB::AdjointAbsVec, adjF::Adjoint{<:Any,<:Factorization}) = adjoint(adjF.parent \ adjB.parent)
 /(B::TransposeAbsVec, adjF::Adjoint{<:Any,<:Factorization}) = adjoint(adjF.parent \ adjoint(B))
 
+
 # support the same 3-arg idiom as in our other in-place A_*_B functions:
 function ldiv!(Y::AbstractVecOrMat, A::Factorization, B::AbstractVecOrMat)
     require_one_based_indexing(Y, B)
diff --git a/stdlib/LinearAlgebra/src/svd.jl b/stdlib/LinearAlgebra/src/svd.jl
index 6991cc9bf13d4..7adf13e56ce58 100644
--- a/stdlib/LinearAlgebra/src/svd.jl
+++ b/stdlib/LinearAlgebra/src/svd.jl
@@ -235,7 +235,7 @@ svdvals(A::AbstractVector{<:BlasFloat}) = [norm(A)]
 svdvals(x::Number) = abs(x)
 svdvals(S::SVD{<:Any,T}) where {T} = (S.S)::Vector{T}
 
-# SVD least squares
+### SVD least squares ###
 function ldiv!(A::SVD{T}, B::StridedVecOrMat) where T
     m, n = size(A)
     k = searchsortedlast(A.S, eps(real(T))*A.S[1], rev=true)

From 63480ba0c1dee73d2eb8ff4afcf9aef60bb3f10c Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreas@noack.dk>
Date: Wed, 26 May 2021 15:09:14 +0200
Subject: [PATCH 11/15] Remove redundant \ definitions in diagonal.jl

---
 stdlib/LinearAlgebra/src/diagonal.jl | 8 --------
 1 file changed, 8 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl
index 204ff56e0c443..464e85facc640 100644
--- a/stdlib/LinearAlgebra/src/diagonal.jl
+++ b/stdlib/LinearAlgebra/src/diagonal.jl
@@ -488,14 +488,6 @@ rdiv!(A::AbstractMatrix{T}, transD::Transpose{<:Any,<:Diagonal{T}}) where {T} =
 (/)(A::Union{StridedMatrix, AbstractTriangular}, D::Diagonal) =
     rdiv!((typeof(oneunit(eltype(D))/oneunit(eltype(A)))).(A), D)
 
-(\)(F::Factorization, D::Diagonal) =
-    ldiv!(F, Matrix{typeof(oneunit(eltype(D))/oneunit(eltype(F)))}(D))
-\(adjF::Adjoint{<:Any,<:Factorization}, D::Diagonal) =
-    (F = adjF.parent; ldiv!(adjoint(F), Matrix{typeof(oneunit(eltype(D))/oneunit(eltype(F)))}(D)))
-(\)(A::Union{QR,QRCompactWY,QRPivoted}, B::Diagonal) =
-    invoke(\, Tuple{Union{QR,QRCompactWY,QRPivoted}, AbstractVecOrMat}, A, B)
-
-
 @inline function kron!(C::AbstractMatrix, A::Diagonal, B::Diagonal)
     valA = A.diag; nA = length(valA)
     valB = B.diag; nB = length(valB)

From 66a3a34b6e19f91b8cea3662e942b11162d0d4cf Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreas@noack.dk>
Date: Thu, 27 May 2021 13:25:26 +0200
Subject: [PATCH 12/15] Add Factorization constructors for SVD

---
 stdlib/LinearAlgebra/src/svd.jl | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/stdlib/LinearAlgebra/src/svd.jl b/stdlib/LinearAlgebra/src/svd.jl
index 7adf13e56ce58..3290b4529be37 100644
--- a/stdlib/LinearAlgebra/src/svd.jl
+++ b/stdlib/LinearAlgebra/src/svd.jl
@@ -72,6 +72,11 @@ function SVD{T}(U::AbstractArray, S::AbstractVector{Tr}, Vt::AbstractArray) wher
         convert(AbstractArray{T}, Vt))
 end
 
+SVD{T}(F::SVD) where {T} = SVD(
+    convert(AbstractMatrix{T}, F.U),
+    convert(AbstractVector{real(T)}, F.S),
+    convert(AbstractMatrix{T}, F.Vt))
+Factorization{T}(F::SVD) where {T} = SVD{T}(F)
 
 # iteration for destructuring into components
 Base.iterate(S::SVD) = (S.U, Val(:S))

From 6af57f345aaedaa000a766ca4f0011907bedf82a Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreas@noack.dk>
Date: Fri, 28 May 2021 08:37:44 +0200
Subject: [PATCH 13/15] Disambiguate between the specialized \ for real
 lhs-complex rhs and then new \ for LAPACKFactorizations.

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

diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl
index 98638fc6de7cf..f94882a1fc259 100644
--- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl
+++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl
@@ -451,6 +451,9 @@ function (\)(F::Union{<:LAPACKFactorizations,Adjoint{<:Any,<:LAPACKFactorization
     # the complete rhs
     return _cut_B(BB, 1:n)
 end
+# disambiguate
+(\)(F::LAPACKFactorizations{T}, B::VecOrMat{Complex{T}}) where {T<:BlasReal} =
+    invoke(\, Tuple{Factorization{T}, VecOrMat{Complex{T}}}, F, B)
 
 """
     LinearAlgebra.peakflops(n::Integer=2000; parallel::Bool=false)

From a6f222c5adb7391815c4e5e07da1539196d96c07 Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreas@noack.dk>
Date: Fri, 28 May 2021 18:39:32 +0200
Subject: [PATCH 14/15] Adjustments based on review

---
 stdlib/LinearAlgebra/src/LinearAlgebra.jl | 3 +--
 stdlib/LinearAlgebra/src/qr.jl            | 2 +-
 2 files changed, 2 insertions(+), 3 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl
index f94882a1fc259..07bb954807361 100644
--- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl
+++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl
@@ -438,7 +438,6 @@ function (\)(F::Union{<:LAPACKFactorizations,Adjoint{<:Any,<:LAPACKFactorization
 
     if n > size(B, 1)
         # Underdetermined
-        fill!(BB, 0)
         copyto!(view(BB, 1:m, :), B)
     else
         copyto!(BB, B)
@@ -446,7 +445,7 @@ function (\)(F::Union{<:LAPACKFactorizations,Adjoint{<:Any,<:LAPACKFactorization
 
     ldiv!(FF, BB)
 
-    # For tall problems, we compute a least sqaures solution so only part
+    # For tall problems, we compute a least squares solution so only part
     # of the rhs should be returned from \ while ldiv! uses (and returns)
     # the complete rhs
     return _cut_B(BB, 1:n)
diff --git a/stdlib/LinearAlgebra/src/qr.jl b/stdlib/LinearAlgebra/src/qr.jl
index a2f8a49639c9e..d0ec430347193 100644
--- a/stdlib/LinearAlgebra/src/qr.jl
+++ b/stdlib/LinearAlgebra/src/qr.jl
@@ -953,7 +953,7 @@ function ldiv!(Fadj::Adjoint{<:Any,<:Union{QR,QRCompactWY,QRPivoted}}, B::Abstra
     require_one_based_indexing(B)
     m, n = size(Fadj)
 
-    # We don't allow solutions overdetermined systems. It would at least be
+    # We don't allow solutions overdetermined systems
     if m > n
         throw(DimensionMismatch("overdetermined systems are not supported"))
     end

From 285faa8b540ea97581684b5274f9b96f90ff6c90 Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreas@noack.dk>
Date: Sat, 29 May 2021 08:39:15 +0200
Subject: [PATCH 15/15] Fixes for new pivoting syntax

---
 stdlib/LinearAlgebra/test/qr.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/stdlib/LinearAlgebra/test/qr.jl b/stdlib/LinearAlgebra/test/qr.jl
index b713d01ab4e8d..dd5a0db40dd95 100644
--- a/stdlib/LinearAlgebra/test/qr.jl
+++ b/stdlib/LinearAlgebra/test/qr.jl
@@ -392,7 +392,7 @@ end
             end
 
             @testset "QRPivoted" begin
-                F = LinearAlgebra.qr(A, Val(true))
+                F = LinearAlgebra.qr(A, ColumnNorm())
                 x = F'\b
                 @test x ≈ A'\b
                 @test length(size(x)) == length(size(b))
@@ -402,8 +402,8 @@ end
         @test_throws DimensionMismatch("arguments must have the same number of rows") qr(randn(n, n + 1))'\b
         @test_throws DimensionMismatch("overdetermined systems are not supported")    LinearAlgebra.qrfactUnblocked!(randn(n - 2, n))'\b
         @test_throws DimensionMismatch("arguments must have the same number of rows") LinearAlgebra.qrfactUnblocked!(randn(n, n + 1))'\b
-        @test_throws DimensionMismatch("overdetermined systems are not supported")    qr(randn(n - 2, n), Val(true))'\b
-        @test_throws DimensionMismatch("arguments must have the same number of rows") qr(randn(n, n + 1), Val(true))'\b
+        @test_throws DimensionMismatch("overdetermined systems are not supported")    qr(randn(n - 2, n), ColumnNorm())'\b
+        @test_throws DimensionMismatch("arguments must have the same number of rows") qr(randn(n, n + 1), ColumnNorm())'\b
     end
 end