Two performance problems with the Q in QR factorizations

[olof3]

n = 300
B = randn(n, n)
D = Diagonal(ones(n))
F = qr(randn(n, 1))
@time F.Q * B # Works as expected
@time F.Q * D # Slow
@time F.Q * I # Slow
@time F.Q[:, :] # Slow

gives

  0.000505 seconds (4 allocations: 705.734 KiB)
  0.251926 seconds (270.00 k allocations: 448.380 MiB, 10.12% gc time)
  0.246129 seconds (270.00 k allocations: 448.380 MiB, 7.77% gc time)
  0.239383 seconds (270.01 k allocations: 448.380 MiB, 8.09% gc time)

There seems to be two problems

1. Multiplying various matrix types with QRCompactQ.
   The same issue was reported for multiplication with UniformScaling in Unexpected performance loss #758.
   I guess the solution is to define the proper methods. Perhaps this would also fix Missing element promotion in I(n)*F.Q' #801.
2. getindex for QRCompactQ
   This problem could be solved by defining the getindex method properly.
   A bonus would be that one could use F.Q[:,:] to retrieve the full Q, rather than the cumbersome F.Q*Matrix(I,m,m) as suggested in the docstring.

For problem 2., the following method seems to work for me, but I guess it needs further consideration and testing.

function getindex(Q::Union{AbstractQ,Adjoint{<:Any,<:AbstractQ}}, ind1, ind2)
    if ind2 isa Colon
        ind2 = 1:size(Q,2)
    end

    X = zeros(eltype(Q), size(Q,1), length(ind2))
    for (j,i) in enumerate(ind2)
        X[i,j] = 1
    end

    lmul!(Q, X)
    
    if ind2 isa Integer
        return X[ind1, 1]
    else
        return X[ind1, :]
    end
end

Perhaps it would also be good to have

getindex(Q::Union{AbstractQ,Adjoint{<:Any,<:AbstractQ}}, ::Colon, ::Colon) = lmul!(Q, Matrix{eltype(Q)}(I, size(Q)))

If the above methods are introduced it seems like the existing method,

function getindex(Q::AbstractQ, i::Integer, j::Integer)
    x = zeros(eltype(Q), size(Q, 1))
    x[i] = 1
    y = zeros(eltype(Q), size(Q, 2))
    y[j] = 1
    return dot(x, lmul!(Q, y))
end

could be removed.

[dkarrasch]

For multiplication, one can perhaps start from the strided matrix version:

# this needs to be special-cased for all specific types, to avoid method ambiguities :-(
(*)(A::AbstractQ, B::SpecialMatrix) = _qlmul(A, B)

function _qlmul(A::AbstractQ, B)
    TAB = promote_type(eltype(A), eltype(B))
    Anew = convert(AbstractMatrix{TAB}, A)
    if size(A.factors, 1) == size(B, 1)
        Bnew = Matrix{TAB}(B)
    elseif size(A.factors, 2) == size(B, 1)
        Bnew = [Matrix{TAB}(B); zeros(TAB, size(A.factors, 1) - size(B,1), size(B, 2))]
    else
        throw(DimensionMismatch("first dimension of matrix must have size either $(size(A.factors, 1)) or $(size(A.factors, 2))"))
    end
    lmul!(Anew, Bnew)
end

The competing method is, for instance, (*)(A::AbstractMatrix, D::Diagonal), so just having (*)(A::AbstractQ, B::AbstractMatrix) will not work out of the box. Furthermore, we'll need the same methods for multiplying with Q from the right.

[ViralBShah]

Only F.Q * D apears to be slow now, but it now takes 36 secs for me!

[dkarrasch]

Uff, that may be the reason for the significantly longer runtimes of the linalg tests. I'll take care of it.