From 86f656d414463ed20b28aa6138dc3e3917b77a84 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Sun, 31 Oct 2021 13:23:06 +0100 Subject: [PATCH] Allow vectors in reflectorApply! (#42874) Co-authored-by: Sheehan Olver --- stdlib/LinearAlgebra/src/generic.jl | 4 ++-- stdlib/LinearAlgebra/src/qr.jl | 6 +++--- stdlib/LinearAlgebra/test/generic.jl | 12 ++++++++++++ 3 files changed, 17 insertions(+), 5 deletions(-) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index 7ab9b6412cb5b..f519cea042c1f 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -1516,9 +1516,9 @@ end end # apply reflector from left -@inline function reflectorApply!(x::AbstractVector, τ::Number, A::AbstractMatrix) +@inline function reflectorApply!(x::AbstractVector, τ::Number, A::AbstractVecOrMat) require_one_based_indexing(x) - m, n = size(A) + m, n = size(A, 1), size(A, 2) if length(x) != m throw(DimensionMismatch("reflector has length $(length(x)), which must match the first dimension of matrix A, $m")) end diff --git a/stdlib/LinearAlgebra/src/qr.jl b/stdlib/LinearAlgebra/src/qr.jl index 31331c9e159f9..671abc00a6cfe 100644 --- a/stdlib/LinearAlgebra/src/qr.jl +++ b/stdlib/LinearAlgebra/src/qr.jl @@ -270,13 +270,13 @@ function qrfactPivotedUnblocked!(A::AbstractMatrix) # Compute reflector of columns j x = view(A, j:m, j) - τj = LinearAlgebra.reflector!(x) + τj = reflector!(x) τ[j] = τj # Update trailing submatrix with reflector - LinearAlgebra.reflectorApply!(x, τj, view(A, j:m, j+1:n)) + reflectorApply!(x, τj, view(A, j:m, j+1:n)) end - return LinearAlgebra.QRPivoted{eltype(A), typeof(A)}(A, τ, piv) + return QRPivoted{eltype(A), typeof(A)}(A, τ, piv) end # LAPACK version diff --git a/stdlib/LinearAlgebra/test/generic.jl b/stdlib/LinearAlgebra/test/generic.jl index 489b96be56019..23edfac4280b7 100644 --- a/stdlib/LinearAlgebra/test/generic.jl +++ b/stdlib/LinearAlgebra/test/generic.jl @@ -256,6 +256,18 @@ end @test_throws DimensionMismatch reflect!([x; x], y, c, s) end +@testset "LinearAlgebra.reflectorApply!" begin + for T in (Float64, ComplexF64) + x = rand(T, 6) + τ = rand(T) + A = rand(T, 6) + B = LinearAlgebra.reflectorApply!(x, τ, copy(A)) + C = LinearAlgebra.reflectorApply!(x, τ, reshape(copy(A), (length(A), 1))) + @test B[1] ≈ C[1] ≈ A[1] - conj(τ)*(A[1] + dot(x[2:end], A[2:end])) + @test B[2:end] ≈ C[2:end] ≈ A[2:end] - conj(τ)*(A[1] + dot(x[2:end], A[2:end]))*x[2:end] + end +end + @testset "LinearAlgebra.axp(b)y! for element type without commutative multiplication" begin α = [1 2; 3 4] β = [5 6; 7 8]