@less mul!(similar(A), A, Q, true, false) leads you to BLAS.symm!('R', Q.uplo, true, Q.data, A, false, similar(A)). It doesn't look like that has a version which accepts transposed A.

Since A' * Q ≈ (Q' * A)', returning that would be one option. How surprising would it be for * to return an Adjoint matrix?

Since A' * Q ≈ (Q' * A)', returning that would be one option. How surprising would it be for * to return an Adjoint matrix?

Not, sure, a bit surprising, but also nice with some more performance. I'm primarily worried about avoiding generic matmul, though.

Here is another example with common types where generic matmul strikes again.

julia> n = 500; A = randn(n,n); B = rand(Bool, n, n); @btime $A*$B;
  301.147 ms (8 allocations: 1.91 MiB)

Some eltype mismatches are copied to avoid this, e.g. there's a method (*)(A::StridedMatrix{<:BlasReal}, B::StridedMatrix{<:BlasFloat}) which helps the first of these, maybe you could extend it to capture Bool without producing ambiguities?

julia> A = rand(Float32,100,100); B = rand(100,100);

julia> @btime $A * $B;  # special method which copies
  36.000 μs (4 allocations: 156.41 KiB)

julia> @btime Float64.($A) * $B;  # copy by hand, same speed & same alloc.
  34.875 μs (4 allocations: 156.41 KiB)

julia> A = rand(Bool,100,100); B = rand(100,100);

julia> @btime $A * $B;  # generic matmul
  490.875 μs (8 allocations: 78.53 KiB)

julia> @btime Float64.($A) * $B;
  38.500 μs (4 allocations: 156.41 KiB)