using LinearAlgebra

function _wide_qr_ldiv!(A::QRPivoted{Complex{T}, Matrix{Complex{T}}, Vector{Complex{T}}, Vector{Int64}}, B::Matrix{Complex{T}}) where T
    m, n = size(A)
    minmn = min(m,n)
    mB, nB = size(B)
    lmul!(adjoint(A.Q), view(B, 1:m, :))
    R = A.R # makes a copy, used as a buffer below
    @inbounds begin
        if n > m # minimum norm solution
            τ = zeros(T,m)
            for k = m:-1:1 # Trapezoid to triangular by elementary operation
                x = view(R, k, [k; m + 1:n])
                τk = LinearAlgebra.reflector!(x)
                τ[k] = conj(τk)
                for i = 1:k - 1
                    vRi = R[i,k]
                    for j = m + 1:n
                        vRi += R[i,j]*x[j - m + 1]'
                    end
                    vRi *= τk
                    R[i,k] -= vRi
                    for j = m + 1:n
                        R[i,j] -= vRi*x[j - m + 1]
                    end
                end
            end
        end
        ldiv!(UpperTriangular(view(R, :, 1:minmn)), view(B, 1:minmn, :))
        if n > m # Apply elementary transformation to solution
            B[m + 1:mB,1:nB] .= zero(T)
            for j = 1:nB
                for k = 1:m
                    vBj = B[k,j]'
                    for i = m + 1:n
                        vBj += B[i,j]'*R[k,i]'
                    end
                    vBj *= τ[k]
                    B[k,j] -= vBj'
                    for i = m + 1:n
                        B[i,j] -= R[k,i]'*vBj'
                    end
                end
            end
        end
    end
    B[A.p,:] = B[1:n,:]
    return B
end

println("\nOriginal Test Case\n")
A = Complex{BigFloat}[1+im 1-im]
B = Complex{BigFloat}[3+im 0]
BI = copy(B) ;
AF = Complex{Float64}[1+im 1-im]
BF = Complex{Float64}[3+im]
BFI = Complex{Float64}[3+im; 0]
xe = AF\BF
x = _wide_qr_ldiv!(qr(A,ColumnNorm()),reshape(BI,2,1))
display(norm(xe-x))
xf = _wide_qr_ldiv!(qr(AF,ColumnNorm()),reshape(BFI,2,1))
display(norm(xe-xf))

println("\n2x1 Random Test Case\n")
c11 = 5*complex(randn(),randn())
c12 = 5*complex(randn(),randn())
c13 = 5*complex(randn(),randn())
c21 = 5*complex(randn(),randn())
c22 = 5*complex(randn(),randn())
c23 = 5*complex(randn(),randn())
d1 = 5*complex(randn(),randn())
d2 = 5*complex(randn(),randn())
C = Complex{BigFloat}[c11 c12; c21 c22]
D = Complex{BigFloat}[d1; 0]
DI = copy(D)
CF = Complex{Float64}[c11 c12; c21 c22]
DF = Complex{Float64}[d1; 0]
DFI = Complex{Float64}[d1; 0]
xe = CF\DF
x = _wide_qr_ldiv!(qr(C,ColumnNorm()),reshape(DI,2,1))
display(norm(xe-x))
xf = _wide_qr_ldiv!(qr(CF,ColumnNorm()),reshape(DFI,2,1))
display(norm(xe-xf))

println("\n3x1 Random Test Case\n")
C = Complex{BigFloat}[c11 c12 c13; c21 c22 c23]
D = Complex{BigFloat}[d1; 0; 0]
DI = copy(D)
CF = Complex{Float64}[c11 c12 c13; c21 c22 c23]
DF = Complex{Float64}[d1; 0]
DFI = Complex{Float64}[d1; 0; 0]
xe = CF\DF
x = _wide_qr_ldiv!(qr(C,ColumnNorm()),reshape(DI,3,1))
display(norm(xe-x))
xf = _wide_qr_ldiv!(qr(CF,ColumnNorm()),reshape(DFI,3,1))
display(norm(xe-xf))

println("\n3x2 Random Test Case\n")
C = Complex{BigFloat}[c11 c12 c13; c21 c22 c23]
D = Complex{BigFloat}[d1 d2; 0 0; 0 0]
DI = copy(D)
CF = Complex{Float64}[c11 c12 c13; c21 c22 c23]
DF = Complex{Float64}[d1 d2; 0 0]
DFI = Complex{Float64}[d1 d2; 0 0; 0 0]
xe = CF\DF
x = _wide_qr_ldiv!(qr(C,ColumnNorm()),reshape(DI,3,2))
display(norm(xe-x))
xf = _wide_qr_ldiv!(qr(CF,ColumnNorm()),reshape(DFI,3,2))
display(norm(xe-xf))