diff --git a/NEWS.md b/NEWS.md index ff1501104e6be..976c88ee65cb8 100644 --- a/NEWS.md +++ b/NEWS.md @@ -51,6 +51,15 @@ Standard library changes #### LinearAlgebra +* `AbstractQ` no longer subtypes to `AbstractMatrix`. Moreover, `adjoint(Q::AbstractQ)` + no longer wraps `Q` in an `Adjoint` type, but instead in an `AdjointQ`, that itself + subtypes `AbstractQ`. This change accounts for the fact that typically `AbstractQ` + instances behave like function-based, matrix-backed linear operators, and hence don't + allow for efficient indexing. Also, many `AbstractQ` types can act on vectors/matrices + of different size, acting like a matrix with context-dependent size. With this change, + `AbstractQ` has a well-defined API that is described in detail in the + [Julia documentation](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#man-linalg-abstractq) + ([#46196]). * Adjoints and transposes of `Factorization` objects are no longer wrapped in `Adjoint` and `Transpose` wrappers, respectively. Instead, they are wrapped in `AdjointFactorization` and `TranposeFactorization` types, which themselves subtype diff --git a/stdlib/LinearAlgebra/docs/src/index.md b/stdlib/LinearAlgebra/docs/src/index.md index 5b5c43da9875d..cd317b4e36df6 100644 --- a/stdlib/LinearAlgebra/docs/src/index.md +++ b/stdlib/LinearAlgebra/docs/src/index.md @@ -150,7 +150,10 @@ julia> sB\x -1.1086956521739126 -1.4565217391304346 ``` -The `\` operation here performs the linear solution. The left-division operator is pretty powerful and it's easy to write compact, readable code that is flexible enough to solve all sorts of systems of linear equations. + +The `\` operation here performs the linear solution. The left-division operator is pretty +powerful and it's easy to write compact, readable code that is flexible enough to solve all +sorts of systems of linear equations. ## Special matrices @@ -309,6 +312,94 @@ Adjoints and transposes of [`Factorization`](@ref) objects are lazily wrapped in `AdjointFactorization` and `TransposeFactorization` objects, respectively. Generically, transpose of real `Factorization`s are wrapped as `AdjointFactorization`. +## [Orthogonal matrices (`AbstractQ`)](@id man-linalg-abstractq) + +Some matrix factorizations generate orthogonal/unitary "matrix" factors. These +factorizations include QR-related factorizations obtained from calls to [`qr`](@ref), i.e., +`QR`, `QRCompactWY` and `QRPivoted`, the Hessenberg factorization obtained from calls to +[`hessenberg`](@ref), and the LQ factorization obtained from [`lq`](@ref). While these +orthogonal/unitary factors admit a matrix representation, their internal representation +is, for performance and memory reasons, different. Hence, they should be rather viewed as +matrix-backed, function-based linear operators. In particular, reading, for instance, a +column of its matrix representation requires running "matrix"-vector multiplication code, +rather than simply reading out data from memory (possibly filling parts of the vector with +structural zeros). Another clear distinction from other, non-triangular matrix types is +that the underlying multiplication code allows for in-place modification during multiplication. +Furthermore, objects of specific `AbstractQ` subtypes as those created via [`qr`](@ref), +[`hessenberg`](@ref) and [`lq`](@ref) can behave like a square or a rectangular matrix +depending on context: + +```julia +julia> using LinearAlgebra + +julia> Q = qr(rand(3,2)).Q +3×3 LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}, Matrix{Float64}} + +julia> Matrix(Q) +3×2 Matrix{Float64}: + -0.320597 0.865734 + -0.765834 -0.475694 + -0.557419 0.155628 + +julia> Q*I +3×3 Matrix{Float64}: + -0.320597 0.865734 -0.384346 + -0.765834 -0.475694 -0.432683 + -0.557419 0.155628 0.815514 + +julia> Q*ones(2) +3-element Vector{Float64}: + 0.5451367118802273 + -1.241527373086654 + -0.40179067589600226 + +julia> Q*ones(3) +3-element Vector{Float64}: + 0.16079054743832022 + -1.674209978965636 + 0.41372375588835797 + +julia> ones(1,2) * Q' +1×3 Matrix{Float64}: + 0.545137 -1.24153 -0.401791 + +julia> ones(1,3) * Q' +1×3 Matrix{Float64}: + 0.160791 -1.67421 0.413724 +``` + +Due to this distinction from dense or structured matrices, the abstract `AbstractQ` type +does not subtype `AbstractMatrix`, but instead has its own type hierarchy. Custom types +that subtype `AbstractQ` can rely on generic fallbacks if the following interface is satisfied. +For example, for + +```julia +struct MyQ{T} <: LinearAlgebra.AbstractQ{T} + # required fields +end +``` + +provide overloads for + +```julia +Base.size(Q::MyQ) # size of corresponding square matrix representation +Base.convert(::Type{AbstractQ{T}}, Q::MyQ) # eltype promotion [optional] +LinearAlgebra.lmul!(Q::MyQ, x::AbstractVecOrMat) # left-multiplication +LinearAlgebra.rmul!(A::AbstractMatrix, Q::MyQ) # right-multiplication +``` + +If `eltype` promotion is not of interest, the `convert` method is unnecessary, since by +default `convert(::Type{AbstractQ{T}}, Q::AbstractQ{T})` returns `Q` itself. +Adjoints of `AbstractQ`-typed objects are lazily wrapped in an `AdjointQ` wrapper type, +which requires its own `LinearAlgebra.lmul!` and `LinearAlgebra.rmul!` methods. Given this +set of methods, any `Q::MyQ` can be used like a matrix, preferably in a multiplicative +context: multiplication via `*` with scalars, vectors and matrices from left and right, +obtaining a matrix representation of `Q` via `Matrix(Q)` (or `Q*I`) and indexing into the +matrix representation all work. In contrast, addition and subtraction as well as more +generally broadcasting over elements in the matrix representation fail because that would +be highly inefficient. For such use cases, consider computing the matrix representation +up front and cache it for future reuse. + ## Standard functions Linear algebra functions in Julia are largely implemented by calling functions from [LAPACK](http://www.netlib.org/lapack/). @@ -505,6 +596,7 @@ four methods defined, for [`Float32`](@ref), [`Float64`](@ref), [`ComplexF32`](@ and [`ComplexF64`](@ref Complex) arrays. ### [BLAS character arguments](@id stdlib-blas-chars) + Many BLAS functions accept arguments that determine whether to transpose an argument (`trans`), which triangle of a matrix to reference (`uplo` or `ul`), whether the diagonal of a triangular matrix can be assumed to @@ -512,18 +604,21 @@ be all ones (`dA`) or which side of a matrix multiplication the input argument belongs on (`side`). The possibilities are: #### [Multiplication order](@id stdlib-blas-side) + | `side` | Meaning | |:-------|:--------------------------------------------------------------------| | `'L'` | The argument goes on the *left* side of a matrix-matrix operation. | | `'R'` | The argument goes on the *right* side of a matrix-matrix operation. | #### [Triangle referencing](@id stdlib-blas-uplo) + | `uplo`/`ul` | Meaning | |:------------|:------------------------------------------------------| | `'U'` | Only the *upper* triangle of the matrix will be used. | | `'L'` | Only the *lower* triangle of the matrix will be used. | #### [Transposition operation](@id stdlib-blas-trans) + | `trans`/`tX` | Meaning | |:-------------|:--------------------------------------------------------| | `'N'` | The input matrix `X` is not transposed or conjugated. | @@ -531,12 +626,12 @@ the input argument belongs on (`side`). The possibilities are: | `'C'` | The input matrix `X` will be conjugated and transposed. | #### [Unit diagonal](@id stdlib-blas-diag) + | `diag`/`dX` | Meaning | |:------------|:----------------------------------------------------------| | `'N'` | The diagonal values of the matrix `X` will be read. | | `'U'` | The diagonal of the matrix `X` is assumed to be all ones. | - ```@docs LinearAlgebra.BLAS LinearAlgebra.BLAS.set_num_threads @@ -582,6 +677,7 @@ and define matrix-vector operations. [Dongarra-1988]: https://dl.acm.org/doi/10.1145/42288.42291 **return a vector** + ```@docs LinearAlgebra.BLAS.gemv! LinearAlgebra.BLAS.gemv(::Any, ::Any, ::Any, ::Any) @@ -611,6 +707,7 @@ LinearAlgebra.BLAS.trsv ``` **return a matrix** + ```@docs LinearAlgebra.BLAS.ger! # xGERU diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index ed5e8d24ba80f..a29c259dae607 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -15,7 +15,7 @@ import Base: USE_BLAS64, abs, acos, acosh, acot, acoth, acsc, acsch, adjoint, as length, log, map, ndims, one, oneunit, parent, permutedims, power_by_squaring, print_matrix, promote_rule, real, round, sec, sech, setindex!, show, similar, sin, sincos, sinh, size, sqrt, strides, stride, tan, tanh, transpose, trunc, typed_hcat, - vec, zero + vec, view, zero using Base: IndexLinear, promote_eltype, promote_op, promote_typeof, @propagate_inbounds, reduce, typed_hvcat, typed_vcat, require_one_based_indexing, splat @@ -431,8 +431,6 @@ include("tridiag.jl") include("triangular.jl") include("factorization.jl") -include("qr.jl") -include("lq.jl") include("eigen.jl") include("svd.jl") include("symmetric.jl") @@ -443,7 +441,10 @@ include("diagonal.jl") include("symmetriceigen.jl") include("bidiag.jl") include("uniformscaling.jl") +include("qr.jl") +include("lq.jl") include("hessenberg.jl") +include("abstractq.jl") include("givens.jl") include("special.jl") include("bitarray.jl") diff --git a/stdlib/LinearAlgebra/src/abstractq.jl b/stdlib/LinearAlgebra/src/abstractq.jl new file mode 100644 index 0000000000000..e591ca69fe429 --- /dev/null +++ b/stdlib/LinearAlgebra/src/abstractq.jl @@ -0,0 +1,628 @@ +# This file is a part of Julia. License is MIT: https://julialang.org/license + +abstract type AbstractQ{T} end + +struct AdjointQ{T,S<:AbstractQ{T}} <: AbstractQ{T} + Q::S +end + +parent(adjQ::AdjointQ) = adjQ.Q +eltype(::Type{<:AbstractQ{T}}) where {T} = T +ndims(::AbstractQ) = 2 + +# inversion/adjoint/transpose +inv(Q::AbstractQ) = Q' +adjoint(Q::AbstractQ) = AdjointQ(Q) +transpose(Q::AbstractQ{<:Real}) = AdjointQ(Q) +transpose(Q::AbstractQ) = error("transpose not implemented for $(typeof(Q)). Consider using adjoint instead of transpose.") +adjoint(adjQ::AdjointQ) = adjQ.Q + +# promotion with AbstractMatrix, at least for equal eltypes +promote_rule(::Type{<:AbstractMatrix{T}}, ::Type{<:AbstractQ{T}}) where {T} = + (@inline; Union{AbstractMatrix{T},AbstractQ{T}}) + +# conversion +AbstractQ{S}(Q::AbstractQ{S}) where {S} = Q +# the following eltype promotion needs to be defined for each subtype +# convert(::Type{AbstractQ{T}}, Q::QType) where {T} = QType{T}(Q) +convert(::Type{AbstractQ{T}}, Q::AbstractQ{T}) where {T} = Q +convert(::Type{AbstractQ{T}}, adjQ::AdjointQ{T}) where {T} = adjQ +convert(::Type{AbstractQ{T}}, adjQ::AdjointQ) where {T} = convert(AbstractQ{T}, adjQ.Q)' + +# ... to matrix +Matrix{T}(Q::AbstractQ) where {T} = convert(Matrix{T}, Q*I) # generic fallback, yields square matrix +Matrix{T}(adjQ::AdjointQ{S}) where {T,S} = convert(Matrix{T}, lmul!(adjQ, Matrix{S}(I, size(adjQ)))) +Matrix(Q::AbstractQ{T}) where {T} = Matrix{T}(Q) +Array{T}(Q::AbstractQ) where {T} = Matrix{T}(Q) +Array(Q::AbstractQ) = Matrix(Q) +convert(::Type{T}, Q::AbstractQ) where {T<:Array} = T(Q) +convert(::Type{T}, Q::AbstractQ) where {T<:Matrix} = T(Q) +# legacy +@deprecate(convert(::Type{AbstractMatrix{T}}, Q::AbstractQ) where {T}, + convert(LinearAlgebra.AbstractQ{T}, Q)) + +function size(Q::AbstractQ, dim::Integer) + if dim < 1 + throw(BoundsError()) + elseif dim <= 2 # && 1 <= dim + return size(Q)[dim] + else # 2 < dim + return 1 + end +end +size(adjQ::AdjointQ) = reverse(size(adjQ.Q)) + +# pseudo-array behaviour, required for indexing with `begin` or `end` +axes(Q::AbstractQ) = map(Base.oneto, size(Q)) +axes(Q::AbstractQ, d::Integer) = d in (1, 2) ? axes(Q)[d] : Base.OneTo(1) + +copymutable(Q::AbstractQ{T}) where {T} = lmul!(Q, Matrix{T}(I, size(Q))) +copy(Q::AbstractQ) = copymutable(Q) + +# getindex +@inline function getindex(Q::AbstractQ, inds...) + @boundscheck Base.checkbounds_indices(Bool, axes(Q), inds) || Base.throw_boundserror(Q, inds) + return _getindex(Q, inds...) +end +@inline getindex(Q::AbstractQ, ::Colon) = copymutable(Q)[:] +@inline getindex(Q::AbstractQ, ::Colon, ::Colon) = copy(Q) + +@inline _getindex(Q::AbstractQ, inds...) = @inbounds copymutable(Q)[inds...] +@inline function _getindex(Q::AbstractQ, ::Colon, J::AbstractVector{<:Integer}) + Y = zeros(eltype(Q), size(Q, 2), length(J)) + @inbounds for (i,j) in enumerate(J) + Y[j,i] = oneunit(eltype(Q)) + end + lmul!(Q, Y) +end +@inline _getindex(Q::AbstractQ, I::AbstractVector{Int}, J::AbstractVector{Int}) = @inbounds Q[:,J][I,:] +@inline function _getindex(Q::AbstractQ, ::Colon, j::Int) + y = zeros(eltype(Q), size(Q, 2)) + y[j] = oneunit(eltype(Q)) + lmul!(Q, y) +end +@inline _getindex(Q::AbstractQ, i::Int, j::Int) = @inbounds Q[:,j][i] + +# needed because AbstractQ does not subtype AbstractMatrix +qr(Q::AbstractQ{T}, arg...; kwargs...) where {T} = qr!(Matrix{_qreltype(T)}(Q), arg...; kwargs...) +lq(Q::AbstractQ{T}, arg...; kwargs...) where {T} = lq!(Matrix{lq_eltype(T)}(Q), arg...; kwargs...) +hessenberg(Q::AbstractQ{T}) where {T} = hessenberg!(Matrix{eigtype(T)}(Q)) + +# needed when used interchangeably with AbstractMatrix (analogous to views of ranges) +view(A::AbstractQ, I...) = getindex(A, I...) + +# specialization avoiding the fallback using slow `getindex` +function copyto!(dest::AbstractMatrix, src::AbstractQ) + copyto!(dest, I) + lmul!(src, dest) +end +# needed to resolve method ambiguities +function copyto!(dest::PermutedDimsArray{T,2,perm}, src::AbstractQ) where {T,perm} + if perm == (1, 2) + copyto!(parent(dest), src) + else + @assert perm == (2, 1) # there are no other permutations of two indices + if T <: Real + copyto!(parent(dest), I) + lmul!(src', parent(dest)) + else + # LAPACK does not offer inplace lmul!(transpose(Q), B) for complex Q + tmp = similar(parent(dest)) + copyto!(tmp, I) + rmul!(tmp, src) + permutedims!(parent(dest), tmp, (2, 1)) + end + end + return dest +end + +function show(io::IO, ::MIME{Symbol("text/plain")}, Q::AbstractQ) + print(io, Base.dims2string(size(Q)), ' ', summary(Q)) +end + +# multiplication +(*)(Q::AbstractQ, J::UniformScaling) = Q*J.λ +function (*)(Q::AbstractQ, b::Number) + T = promote_type(eltype(Q), typeof(b)) + lmul!(convert(AbstractQ{T}, Q), Matrix{T}(b*I, size(Q))) +end +function (*)(A::AbstractQ, B::AbstractVecOrMat) + T = promote_type(eltype(A), eltype(B)) + lmul!(convert(AbstractQ{T}, A), copy_similar(B, T)) +end + +(*)(J::UniformScaling, Q::AbstractQ) = J.λ*Q +function (*)(a::Number, Q::AbstractQ) + T = promote_type(typeof(a), eltype(Q)) + rmul!(Matrix{T}(a*I, size(Q)), convert(AbstractQ{T}, Q)) +end +*(a::AbstractVector, Q::AbstractQ) = reshape(a, length(a), 1) * Q +function (*)(A::AbstractMatrix, Q::AbstractQ) + T = promote_type(eltype(A), eltype(Q)) + return rmul!(copy_similar(A, T), convert(AbstractQ{T}, Q)) +end +(*)(u::AdjointAbsVec, Q::AbstractQ) = (Q'u')' + +# AbstractQ * Triangular +lmul!(Q::AbstractQ, B::AbstractTriangular) = lmul!(Q, full!(B)) +rmul!(A::AbstractTriangular, Q::AbstractQ) = rmul!(full!(A), Q) + +### Q*Q (including adjoints) +*(Q::AbstractQ, P::AbstractQ) = Q * (P*I) + +### mul! +function mul!(C::AbstractVecOrMat{T}, Q::AbstractQ{T}, B::Union{AbstractVecOrMat{T},AbstractQ{T}}) where {T} + require_one_based_indexing(C, B) + mB = size(B, 1) + mC = size(C, 1) + if mB < mC + inds = CartesianIndices(axes(B)) + copyto!(view(C, inds), B) + C[CartesianIndices((mB+1:mC, axes(C, 2)))] .= zero(T) + return lmul!(Q, C) + else + return lmul!(Q, copyto!(C, B)) + end +end +mul!(C::AbstractVecOrMat{T}, A::AbstractVecOrMat{T}, Q::AbstractQ{T}) where {T} = rmul!(copyto!(C, A), Q) +mul!(C::AbstractVecOrMat{T}, adjQ::AdjointQ{T}, B::AbstractVecOrMat{T}) where {T} = lmul!(adjQ, copyto!(C, B)) +mul!(C::AbstractVecOrMat{T}, A::AbstractVecOrMat{T}, adjQ::AdjointQ{T}) where {T} = rmul!(copyto!(C, A), adjQ) + +### division +\(Q::AbstractQ, A::AbstractVecOrMat) = Q'*A +/(A::AbstractVecOrMat, Q::AbstractQ) = A*Q' +ldiv!(Q::AbstractQ, A::AbstractVecOrMat) = lmul!(Q', A) +ldiv!(C::AbstractVecOrMat, Q::AbstractQ, A::AbstractVecOrMat) = mul!(C, Q', A) +rdiv!(A::AbstractVecOrMat, Q::AbstractQ) = rmul!(A, Q') + +logabsdet(Q::AbstractQ) = (d = det(Q); return log(abs(d)), sign(d)) +function logdet(A::AbstractQ) + d, s = logabsdet(A) + return d + log(s) +end + +########################################################### +################ Q from QR decompositions ################# +########################################################### + +""" + QRPackedQ <: LinearAlgebra.AbstractQ + +The orthogonal/unitary ``Q`` matrix of a QR factorization stored in [`QR`](@ref) or +[`QRPivoted`](@ref) format. +""" +struct QRPackedQ{T,S<:AbstractMatrix{T},C<:AbstractVector{T}} <: AbstractQ{T} + factors::S + τ::C + + function QRPackedQ{T,S,C}(factors, τ) where {T,S<:AbstractMatrix{T},C<:AbstractVector{T}} + require_one_based_indexing(factors, τ) + new{T,S,C}(factors, τ) + end +end +QRPackedQ(factors::AbstractMatrix{T}, τ::AbstractVector{T}) where {T} = + QRPackedQ{T,typeof(factors),typeof(τ)}(factors, τ) +QRPackedQ{T}(factors::AbstractMatrix, τ::AbstractVector) where {T} = + QRPackedQ(convert(AbstractMatrix{T}, factors), convert(AbstractVector{T}, τ)) +# backwards-compatible constructors (remove with Julia 2.0) +@deprecate(QRPackedQ{T,S}(factors::AbstractMatrix{T}, τ::AbstractVector{T}) where {T,S}, + QRPackedQ{T,S,typeof(τ)}(factors, τ), false) + +""" + QRCompactWYQ <: LinearAlgebra.AbstractQ + +The orthogonal/unitary ``Q`` matrix of a QR factorization stored in [`QRCompactWY`](@ref) +format. +""" +struct QRCompactWYQ{S, M<:AbstractMatrix{S}, C<:AbstractMatrix{S}} <: AbstractQ{S} + factors::M + T::C + + function QRCompactWYQ{S,M,C}(factors, T) where {S,M<:AbstractMatrix{S},C<:AbstractMatrix{S}} + require_one_based_indexing(factors, T) + new{S,M,C}(factors, T) + end +end +QRCompactWYQ(factors::AbstractMatrix{S}, T::AbstractMatrix{S}) where {S} = + QRCompactWYQ{S,typeof(factors),typeof(T)}(factors, T) +QRCompactWYQ{S}(factors::AbstractMatrix, T::AbstractMatrix) where {S} = + QRCompactWYQ(convert(AbstractMatrix{S}, factors), convert(AbstractMatrix{S}, T)) +# backwards-compatible constructors (remove with Julia 2.0) +@deprecate(QRCompactWYQ{S,M}(factors::AbstractMatrix{S}, T::AbstractMatrix{S}) where {S,M}, + QRCompactWYQ{S,M,typeof(T)}(factors, T), false) + +QRPackedQ{T}(Q::QRPackedQ) where {T} = QRPackedQ(convert(AbstractMatrix{T}, Q.factors), convert(Vector{T}, Q.τ)) +QRCompactWYQ{S}(Q::QRCompactWYQ) where {S} = QRCompactWYQ(convert(AbstractMatrix{S}, Q.factors), convert(AbstractMatrix{S}, Q.T)) + +AbstractQ{S}(Q::QRPackedQ) where {S} = QRPackedQ{S}(Q) +AbstractQ{S}(Q::QRCompactWYQ) where {S} = QRCompactWYQ{S}(Q) +# override generic square fallback +Matrix{T}(Q::Union{QRCompactWYQ{S},QRPackedQ{S}}) where {T,S} = + convert(Matrix{T}, lmul!(Q, Matrix{S}(I, size(Q, 1), min(size(Q.factors)...)))) +Matrix(Q::Union{QRCompactWYQ{S},QRPackedQ{S}}) where {S} = Matrix{S}(Q) + +convert(::Type{AbstractQ{T}}, Q::QRPackedQ) where {T} = QRPackedQ{T}(Q) +convert(::Type{AbstractQ{T}}, Q::QRCompactWYQ) where {T} = QRCompactWYQ{T}(Q) + +size(Q::Union{QRCompactWYQ,QRPackedQ}, dim::Integer) = + size(Q.factors, dim == 2 ? 1 : dim) +size(Q::Union{QRCompactWYQ,QRPackedQ}) = (n = size(Q.factors, 1); (n, n)) + +## Multiplication +### QB +lmul!(A::QRCompactWYQ{T,<:StridedMatrix}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = + LAPACK.gemqrt!('L', 'N', A.factors, A.T, B) +lmul!(A::QRPackedQ{T,<:StridedMatrix}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = + LAPACK.ormqr!('L', 'N', A.factors, A.τ, B) +function lmul!(A::QRPackedQ, B::AbstractVecOrMat) + require_one_based_indexing(B) + mA, nA = size(A.factors) + mB, nB = size(B,1), size(B,2) + if mA != mB + throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)")) + end + Afactors = A.factors + @inbounds begin + for k = min(mA,nA):-1:1 + for j = 1:nB + vBj = B[k,j] + for i = k+1:mB + vBj += conj(Afactors[i,k])*B[i,j] + end + vBj = A.τ[k]*vBj + B[k,j] -= vBj + for i = k+1:mB + B[i,j] -= Afactors[i,k]*vBj + end + end + end + end + B +end +lmul!(Q::QRPackedQ, B::AbstractTriangular) = lmul!(Q, full!(B)) # disambiguation + +### QcB +lmul!(adjQ::AdjointQ{<:Any,<:QRCompactWYQ{T,<:StridedMatrix}}, B::StridedVecOrMat{T}) where {T<:BlasReal} = + (Q = adjQ.Q; LAPACK.gemqrt!('L', 'T', Q.factors, Q.T, B)) +lmul!(adjQ::AdjointQ{<:Any,<:QRCompactWYQ{T,<:StridedMatrix}}, B::StridedVecOrMat{T}) where {T<:BlasComplex} = + (Q = adjQ.Q; LAPACK.gemqrt!('L', 'C', Q.factors, Q.T, B)) +lmul!(adjQ::AdjointQ{<:Any,<:QRPackedQ{T,<:StridedMatrix}}, B::StridedVecOrMat{T}) where {T<:BlasReal} = + (Q = adjQ.Q; LAPACK.ormqr!('L', 'T', Q.factors, Q.τ, B)) +lmul!(adjQ::AdjointQ{<:Any,<:QRPackedQ{T,<:StridedMatrix}}, B::StridedVecOrMat{T}) where {T<:BlasComplex} = + (Q = adjQ.Q; LAPACK.ormqr!('L', 'C', Q.factors, Q.τ, B)) +function lmul!(adjA::AdjointQ{<:Any,<:QRPackedQ}, B::AbstractVecOrMat) + require_one_based_indexing(B) + A = adjA.Q + mA, nA = size(A.factors) + mB, nB = size(B,1), size(B,2) + if mA != mB + throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)")) + end + Afactors = A.factors + @inbounds begin + for k = 1:min(mA,nA) + for j = 1:nB + vBj = B[k,j] + for i = k+1:mB + vBj += conj(Afactors[i,k])*B[i,j] + end + vBj = conj(A.τ[k])*vBj + B[k,j] -= vBj + for i = k+1:mB + B[i,j] -= Afactors[i,k]*vBj + end + end + end + end + B +end +lmul!(Q::AdjointQ{<:Any,<:QRPackedQ}, B::AbstractTriangular) = lmul!(Q, full!(B)) # disambiguation + +### AQ +rmul!(A::StridedVecOrMat{T}, B::QRCompactWYQ{T,<:StridedMatrix}) where {T<:BlasFloat} = + LAPACK.gemqrt!('R', 'N', B.factors, B.T, A) +rmul!(A::StridedVecOrMat{T}, B::QRPackedQ{T,<:StridedMatrix}) where {T<:BlasFloat} = + LAPACK.ormqr!('R', 'N', B.factors, B.τ, A) +function rmul!(A::AbstractMatrix, Q::QRPackedQ) + require_one_based_indexing(A) + mQ, nQ = size(Q.factors) + mA, nA = size(A,1), size(A,2) + if nA != mQ + throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but matrix Q has dimensions ($mQ, $nQ)")) + end + Qfactors = Q.factors + @inbounds begin + for k = 1:min(mQ,nQ) + for i = 1:mA + vAi = A[i,k] + for j = k+1:mQ + vAi += A[i,j]*Qfactors[j,k] + end + vAi = vAi*Q.τ[k] + A[i,k] -= vAi + for j = k+1:nA + A[i,j] -= vAi*conj(Qfactors[j,k]) + end + end + end + end + A +end +rmul!(A::AbstractTriangular, Q::QRPackedQ) = rmul!(full!(A), Q) # disambiguation + +### AQc +rmul!(A::StridedVecOrMat{T}, adjQ::AdjointQ{<:Any,<:QRCompactWYQ{T}}) where {T<:BlasReal} = + (Q = adjQ.Q; LAPACK.gemqrt!('R', 'T', Q.factors, Q.T, A)) +rmul!(A::StridedVecOrMat{T}, adjQ::AdjointQ{<:Any,<:QRCompactWYQ{T}}) where {T<:BlasComplex} = + (Q = adjQ.Q; LAPACK.gemqrt!('R', 'C', Q.factors, Q.T, A)) +rmul!(A::StridedVecOrMat{T}, adjQ::AdjointQ{<:Any,<:QRPackedQ{T}}) where {T<:BlasReal} = + (Q = adjQ.Q; LAPACK.ormqr!('R', 'T', Q.factors, Q.τ, A)) +rmul!(A::StridedVecOrMat{T}, adjQ::AdjointQ{<:Any,<:QRPackedQ{T}}) where {T<:BlasComplex} = + (Q = adjQ.Q; LAPACK.ormqr!('R', 'C', Q.factors, Q.τ, A)) +function rmul!(A::AbstractMatrix, adjQ::AdjointQ{<:Any,<:QRPackedQ}) + require_one_based_indexing(A) + Q = adjQ.Q + mQ, nQ = size(Q.factors) + mA, nA = size(A,1), size(A,2) + if nA != mQ + throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but matrix Q has dimensions ($mQ, $nQ)")) + end + Qfactors = Q.factors + @inbounds begin + for k = min(mQ,nQ):-1:1 + for i = 1:mA + vAi = A[i,k] + for j = k+1:mQ + vAi += A[i,j]*Qfactors[j,k] + end + vAi = vAi*conj(Q.τ[k]) + A[i,k] -= vAi + for j = k+1:nA + A[i,j] -= vAi*conj(Qfactors[j,k]) + end + end + end + end + A +end +rmul!(A::AbstractTriangular, Q::AdjointQ{<:Any,<:QRPackedQ}) = rmul!(full!(A), Q) # disambiguation + +det(Q::QRPackedQ) = _det_tau(Q.τ) +det(Q::QRCompactWYQ) = + prod(i -> _det_tau(_diagview(Q.T[:, i:min(i + size(Q.T, 1), size(Q.T, 2))])), + 1:size(Q.T, 1):size(Q.T, 2)) + +_diagview(A) = @view A[diagind(A)] + +# Compute `det` from the number of Householder reflections. Handle +# the case `Q.τ` contains zeros. +_det_tau(τs::AbstractVector{<:Real}) = + isodd(count(!iszero, τs)) ? -one(eltype(τs)) : one(eltype(τs)) + +# In complex case, we need to compute the non-unit eigenvalue `λ = 1 - c*τ` +# (where `c = v'v`) of each Householder reflector. As we know that the +# reflector must have the determinant of 1, it must satisfy `abs2(λ) == 1`. +# Combining this with the constraint `c > 0`, it turns out that the eigenvalue +# (hence the determinant) can be computed as `λ = -sign(τ)^2`. +# See: https://github.com/JuliaLang/julia/pull/32887#issuecomment-521935716 +_det_tau(τs) = prod(τ -> iszero(τ) ? one(τ) : -sign(τ)^2, τs) + +########################################################### +######## Q from Hessenberg decomposition ################## +########################################################### + +""" + HessenbergQ <: AbstractQ + +Given a [`Hessenberg`](@ref) factorization object `F`, `F.Q` returns +a `HessenbergQ` object, which is an implicit representation of the unitary +matrix `Q` in the Hessenberg factorization `QHQ'` represented by `F`. +This `F.Q` object can be efficiently multiplied by matrices or vectors, +and can be converted to an ordinary matrix type with `Matrix(F.Q)`. +""" +struct HessenbergQ{T,S<:AbstractMatrix,W<:AbstractVector,sym} <: AbstractQ{T} + uplo::Char + factors::S + τ::W + function HessenbergQ{T,S,W,sym}(uplo::AbstractChar, factors, τ) where {T,S<:AbstractMatrix,W<:AbstractVector,sym} + new(uplo, factors, τ) + end +end +HessenbergQ(F::Hessenberg{<:Any,<:UpperHessenberg,S,W}) where {S,W} = HessenbergQ{eltype(F.factors),S,W,false}(F.uplo, F.factors, F.τ) +HessenbergQ(F::Hessenberg{<:Any,<:SymTridiagonal,S,W}) where {S,W} = HessenbergQ{eltype(F.factors),S,W,true}(F.uplo, F.factors, F.τ) + +size(Q::HessenbergQ, dim::Integer) = size(getfield(Q, :factors), dim == 2 ? 1 : dim) +size(Q::HessenbergQ) = size(Q, 1), size(Q, 2) + +# HessenbergQ from LAPACK/BLAS (as opposed to Julia libraries like GenericLinearAlgebra) +const BlasHessenbergQ{T,sym} = HessenbergQ{T,<:StridedMatrix{T},<:StridedVector{T},sym} where {T<:BlasFloat,sym} + +## reconstruct the original matrix +Matrix{T}(Q::BlasHessenbergQ{<:Any,false}) where {T} = convert(Matrix{T}, LAPACK.orghr!(1, size(Q.factors, 1), copy(Q.factors), Q.τ)) +Matrix{T}(Q::BlasHessenbergQ{<:Any,true}) where {T} = convert(Matrix{T}, LAPACK.orgtr!(Q.uplo, copy(Q.factors), Q.τ)) + +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::StridedVecOrMat{T}, Q::BlasHessenbergQ{T,false}) where {T<:BlasFloat} = + LAPACK.ormhr!('R', 'N', 1, size(Q.factors, 1), Q.factors, Q.τ, X) +lmul!(adjQ::AdjointQ{<:Any,<:BlasHessenbergQ{T,false}}, X::StridedVecOrMat{T}) where {T<:BlasFloat} = + (Q = adjQ.Q; LAPACK.ormhr!('L', ifelse(T<:Real, 'T', 'C'), 1, size(Q.factors, 1), Q.factors, Q.τ, X)) +rmul!(X::StridedVecOrMat{T}, adjQ::AdjointQ{<:Any,<:BlasHessenbergQ{T,false}}) where {T<:BlasFloat} = + (Q = adjQ.Q; 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::StridedVecOrMat{T}, Q::BlasHessenbergQ{T,true}) where {T<:BlasFloat} = + LAPACK.ormtr!('R', Q.uplo, 'N', Q.factors, Q.τ, X) +lmul!(adjQ::AdjointQ{<:Any,<:BlasHessenbergQ{T,true}}, X::StridedVecOrMat{T}) where {T<:BlasFloat} = + (Q = adjQ.Q; LAPACK.ormtr!('L', Q.uplo, ifelse(T<:Real, 'T', 'C'), Q.factors, Q.τ, X)) +rmul!(X::StridedVecOrMat{T}, adjQ::AdjointQ{<:Any,<:BlasHessenbergQ{T,true}}) where {T<:BlasFloat} = + (Q = adjQ.Q; 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,<:StridedVecOrMat{T}}, Q::HessenbergQ{T}) where {T} = lmul!(Q', X')' +lmul!(adjQ::AdjointQ{<:Any,<:HessenbergQ{T}}, X::Adjoint{T,<:StridedVecOrMat{T}}) where {T} = rmul!(X', adjQ')' +rmul!(X::Adjoint{T,<:StridedVecOrMat{T}}, adjQ::AdjointQ{<:Any,<:HessenbergQ{T}}) where {T} = lmul!(adjQ', X')' + +# flexible left-multiplication (and adjoint right-multiplication) +function (*)(Q::Union{QRPackedQ,QRCompactWYQ,HessenbergQ}, b::AbstractVector) + T = promote_type(eltype(Q), eltype(b)) + if size(Q.factors, 1) == length(b) + bnew = copy_similar(b, T) + elseif size(Q.factors, 2) == length(b) + bnew = [b; zeros(T, size(Q.factors, 1) - length(b))] + else + throw(DimensionMismatch("vector must have length either $(size(Q.factors, 1)) or $(size(Q.factors, 2))")) + end + lmul!(convert(AbstractQ{T}, Q), bnew) +end +function (*)(Q::Union{QRPackedQ,QRCompactWYQ,HessenbergQ}, B::AbstractMatrix) + T = promote_type(eltype(Q), eltype(B)) + if size(Q.factors, 1) == size(B, 1) + Bnew = copy_similar(B, T) + elseif size(Q.factors, 2) == size(B, 1) + Bnew = [B; zeros(T, size(Q.factors, 1) - size(B,1), size(B, 2))] + else + throw(DimensionMismatch("first dimension of matrix must have size either $(size(Q.factors, 1)) or $(size(Q.factors, 2))")) + end + lmul!(convert(AbstractQ{T}, Q), Bnew) +end +function (*)(A::AbstractMatrix, adjQ::AdjointQ{<:Any,<:Union{QRPackedQ,QRCompactWYQ,HessenbergQ}}) + Q = adjQ.Q + T = promote_type(eltype(A), eltype(adjQ)) + adjQQ = convert(AbstractQ{T}, adjQ) + if size(A, 2) == size(Q.factors, 1) + AA = copy_similar(A, T) + return rmul!(AA, adjQQ) + elseif size(A, 2) == size(Q.factors, 2) + return rmul!([A zeros(T, size(A, 1), size(Q.factors, 1) - size(Q.factors, 2))], adjQQ) + else + throw(DimensionMismatch("matrix A has dimensions $(size(A)) but Q-matrix B has dimensions $(size(adjQ))")) + end +end +(*)(u::AdjointAbsVec, Q::AdjointQ{<:Any,<:Union{QRPackedQ,QRCompactWYQ,HessenbergQ}}) = (Q'u')' + +det(Q::HessenbergQ) = _det_tau(Q.τ) + +########################################################### +################ Q from LQ decomposition ################## +########################################################### + +struct LQPackedQ{T,S<:AbstractMatrix{T},C<:AbstractVector{T}} <: AbstractQ{T} + factors::S + τ::C +end + +LQPackedQ{T}(Q::LQPackedQ) where {T} = LQPackedQ(convert(AbstractMatrix{T}, Q.factors), convert(Vector{T}, Q.τ)) +@deprecate(AbstractMatrix{T}(Q::LQPackedQ) where {T}, + convert(AbstractQ{T}, Q), + false) +Matrix{T}(A::LQPackedQ) where {T} = convert(Matrix{T}, LAPACK.orglq!(copy(A.factors), A.τ)) +convert(::Type{AbstractQ{T}}, Q::LQPackedQ) where {T} = LQPackedQ{T}(Q) + +# size(Q::LQPackedQ) yields the shape of Q's square form +size(Q::LQPackedQ) = (n = size(Q.factors, 2); return n, n) + +## Multiplication +### QB / QcB +lmul!(A::LQPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = LAPACK.ormlq!('L','N',A.factors,A.τ,B) +lmul!(adjA::AdjointQ{<:Any,<:LQPackedQ{T}}, B::StridedVecOrMat{T}) where {T<:BlasReal} = + (A = adjA.Q; LAPACK.ormlq!('L', 'T', A.factors, A.τ, B)) +lmul!(adjA::AdjointQ{<:Any,<:LQPackedQ{T}}, B::StridedVecOrMat{T}) where {T<:BlasComplex} = + (A = adjA.Q; LAPACK.ormlq!('L', 'C', A.factors, A.τ, B)) + +function (*)(adjA::AdjointQ{<:Any,<:LQPackedQ}, B::AbstractVector) + A = adjA.Q + T = promote_type(eltype(A), eltype(B)) + if length(B) == size(A.factors, 2) + C = copy_similar(B, T) + elseif length(B) == size(A.factors, 1) + C = [B; zeros(T, size(A.factors, 2) - size(A.factors, 1), size(B, 2))] + else + throw(DimensionMismatch("length of B, $(length(B)), must equal one of the dimensions of A, $(size(A))")) + end + lmul!(convert(AbstractQ{T}, adjA), C) +end +function (*)(adjA::AdjointQ{<:Any,<:LQPackedQ}, B::AbstractMatrix) + A = adjA.Q + T = promote_type(eltype(A), eltype(B)) + if size(B,1) == size(A.factors,2) + C = copy_similar(B, T) + elseif size(B,1) == size(A.factors,1) + C = [B; zeros(T, size(A.factors, 2) - size(A.factors, 1), size(B, 2))] + else + throw(DimensionMismatch("first dimension of B, $(size(B,1)), must equal one of the dimensions of A, $(size(A))")) + end + lmul!(convert(AbstractQ{T}, adjA), C) +end + +# in-place right-application of LQPackedQs +# these methods require that the applied-to matrix's (A's) number of columns +# match the number of columns (nQ) of the LQPackedQ (Q) (necessary for in-place +# operation, and the underlying LAPACK routine (ormlq) treats the implicit Q +# as its (nQ-by-nQ) square form) +rmul!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasFloat} = + LAPACK.ormlq!('R', 'N', B.factors, B.τ, A) +rmul!(A::StridedMatrix{T}, adjB::AdjointQ{<:Any,<:LQPackedQ{T}}) where {T<:BlasReal} = + (B = adjB.Q; LAPACK.ormlq!('R', 'T', B.factors, B.τ, A)) +rmul!(A::StridedMatrix{T}, adjB::AdjointQ{<:Any,<:LQPackedQ{T}}) where {T<:BlasComplex} = + (B = adjB.Q; LAPACK.ormlq!('R', 'C', B.factors, B.τ, A)) + +# out-of-place right application of LQPackedQs +# +# these methods: (1) check whether the applied-to matrix's (A's) appropriate dimension +# (columns for A_*, rows for Ac_*) matches the number of columns (nQ) of the LQPackedQ (Q), +# and if so effectively apply Q's square form to A without additional shenanigans; and +# (2) if the preceding dimensions do not match, check whether the appropriate dimension of +# A instead matches the number of rows of the matrix of which Q is a factor (i.e. +# size(Q.factors, 1)), and if so implicitly apply Q's truncated form to A by zero extending +# A as necessary for check (1) to pass (if possible) and then applying Q's square form +# +function (*)(A::AbstractVector, Q::LQPackedQ) + T = promote_type(eltype(A), eltype(Q)) + if 1 == size(Q.factors, 2) + C = copy_similar(A, T) + elseif 1 == size(Q.factors, 1) + C = zeros(T, length(A), size(Q.factors, 2)) + copyto!(C, 1, A, 1, length(A)) + else + _rightappdimmismatch("columns") + end + return rmul!(C, convert(AbstractQ{T}, Q)) +end +function (*)(A::AbstractMatrix, Q::LQPackedQ) + T = promote_type(eltype(A), eltype(Q)) + if size(A, 2) == size(Q.factors, 2) + C = copy_similar(A, T) + elseif size(A, 2) == size(Q.factors, 1) + C = zeros(T, size(A, 1), size(Q.factors, 2)) + copyto!(C, 1, A, 1, length(A)) + else + _rightappdimmismatch("columns") + end + return rmul!(C, convert(AbstractQ{T}, Q)) +end +function (*)(adjA::AdjointAbsMat, Q::LQPackedQ) + A = adjA.parent + T = promote_type(eltype(A), eltype(Q)) + if size(A, 1) == size(Q.factors, 2) + C = copy_similar(adjA, T) + elseif size(A, 1) == size(Q.factors, 1) + C = zeros(T, size(A, 2), size(Q.factors, 2)) + adjoint!(view(C, :, 1:size(A, 1)), A) + else + _rightappdimmismatch("rows") + end + return rmul!(C, convert(AbstractQ{T}, Q)) +end +(*)(u::AdjointAbsVec, Q::LQPackedQ) = (Q'u')' + +_rightappdimmismatch(rowsorcols) = + throw(DimensionMismatch(string("the number of $(rowsorcols) of the matrix on the left ", + "must match either (1) the number of columns of the (LQPackedQ) matrix on the right ", + "or (2) the number of rows of that (LQPackedQ) matrix's internal representation ", + "(the factorization's originating matrix's number of rows)"))) + +# In LQ factorization, `Q` is expressed as the product of the adjoint of the +# reflectors. Thus, `det` has to be conjugated. +det(Q::LQPackedQ) = conj(_det_tau(Q.τ)) diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl index 03c5a2bbdeba4..9bb7219965d5d 100644 --- a/stdlib/LinearAlgebra/src/diagonal.jl +++ b/stdlib/LinearAlgebra/src/diagonal.jl @@ -281,15 +281,11 @@ end rmul!(A::AbstractMatrix, D::Diagonal) = @inline mul!(A, A, D) lmul!(D::Diagonal, B::AbstractVecOrMat) = @inline mul!(B, D, B) -function *(A::AdjOrTransAbsMat, D::Diagonal) +function (*)(A::AdjOrTransAbsMat, D::Diagonal) Ac = copy_similar(A, promote_op(*, eltype(A), eltype(D.diag))) rmul!(Ac, D) end - -*(D::Diagonal, adjQ::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) = - rmul!(Array{promote_type(eltype(D), eltype(adjQ))}(D), adjQ) - -function *(D::Diagonal, A::AdjOrTransAbsMat) +function (*)(D::Diagonal, A::AdjOrTransAbsMat) Ac = copy_similar(A, promote_op(*, eltype(A), eltype(D.diag))) lmul!(D, Ac) end diff --git a/stdlib/LinearAlgebra/src/hessenberg.jl b/stdlib/LinearAlgebra/src/hessenberg.jl index 44cd0873c1ee5..ef7570731197a 100644 --- a/stdlib/LinearAlgebra/src/hessenberg.jl +++ b/stdlib/LinearAlgebra/src/hessenberg.jl @@ -167,29 +167,6 @@ function \(U::UnitUpperTriangular, H::UpperHessenberg) UpperHessenberg(HH) end -function *(H::UpperHessenberg, B::Bidiagonal) - TS = promote_op(matprod, eltype(H), eltype(B)) - A = mul!(similar(H, TS, size(H)), H, B) - return B.uplo == 'U' ? UpperHessenberg(A) : A -end -function *(B::Bidiagonal, H::UpperHessenberg) - TS = promote_op(matprod, eltype(B), eltype(H)) - A = mul!(similar(H, TS, size(H)), B, H) - return B.uplo == 'U' ? UpperHessenberg(A) : A -end - -function /(H::UpperHessenberg, B::Bidiagonal) - T = typeof(oneunit(eltype(H))/oneunit(eltype(B))) - A = _rdiv!(similar(H, T, size(H)), H, B) - return B.uplo == 'U' ? UpperHessenberg(A) : A -end - -function \(B::Bidiagonal, H::UpperHessenberg) - T = typeof(oneunit(eltype(B))\oneunit(eltype(H))) - A = ldiv!(similar(H, T, size(H)), B, H) - return B.uplo == 'U' ? UpperHessenberg(A) : A -end - # Solving (H+µI)x = b: we can do this in O(m²) time and O(m) memory # (in-place in x) by the RQ algorithm from: # @@ -482,10 +459,7 @@ julia> A = [4. 9. 7.; 4. 4. 1.; 4. 3. 2.] julia> F = hessenberg(A) Hessenberg{Float64, UpperHessenberg{Float64, Matrix{Float64}}, Matrix{Float64}, Vector{Float64}, Bool} Q factor: -3×3 LinearAlgebra.HessenbergQ{Float64, Matrix{Float64}, Vector{Float64}, false}: - 1.0 0.0 0.0 - 0.0 -0.707107 -0.707107 - 0.0 -0.707107 0.707107 +3×3 LinearAlgebra.HessenbergQ{Float64, Matrix{Float64}, Vector{Float64}, false} H factor: 3×3 UpperHessenberg{Float64, Matrix{Float64}}: 4.0 -11.3137 -1.41421 @@ -518,43 +492,14 @@ function show(io::IO, mime::MIME"text/plain", F::Hessenberg) show(io, mime, F.H) end -""" - HessenbergQ <: AbstractQ - -Given a [`Hessenberg`](@ref) factorization object `F`, `F.Q` returns -a `HessenbergQ` object, which is an implicit representation of the unitary -matrix `Q` in the Hessenberg factorization `QHQ'` represented by `F`. -This `F.Q` object can be efficiently multiplied by matrices or vectors, -and can be converted to an ordinary matrix type with `Matrix(F.Q)`. -""" -struct HessenbergQ{T,S<:AbstractMatrix,W<:AbstractVector,sym} <: AbstractQ{T} - uplo::Char - factors::S - τ::W - function HessenbergQ{T,S,W,sym}(uplo::AbstractChar, factors, τ) where {T,S<:AbstractMatrix,W<:AbstractVector,sym} - new(uplo, factors, τ) - end -end -HessenbergQ(F::Hessenberg{<:Any,<:UpperHessenberg,S,W}) where {S,W} = HessenbergQ{eltype(F.factors),S,W,false}(F.uplo, F.factors, F.τ) -HessenbergQ(F::Hessenberg{<:Any,<:SymTridiagonal,S,W}) where {S,W} = HessenbergQ{eltype(F.factors),S,W,true}(F.uplo, F.factors, F.τ) - function getproperty(F::Hessenberg, d::Symbol) d === :Q && return HessenbergQ(F) return getfield(F, d) end -size(Q::HessenbergQ, dim::Integer) = size(getfield(Q, :factors), dim == 2 ? 1 : dim) -size(Q::HessenbergQ) = size(Q, 1), size(Q, 2) - Base.propertynames(F::Hessenberg, private::Bool=false) = (:Q, :H, :μ, (private ? (:τ, :factors, :uplo) : ())...) -# HessenbergQ from LAPACK/BLAS (as opposed to Julia libraries like GenericLinearAlgebra) -const BlasHessenbergQ{T,sym} = HessenbergQ{T,<:StridedMatrix{T},<:StridedVector{T},sym} where {T<:BlasFloat,sym} - -## reconstruct the original matrix -Matrix{T}(Q::BlasHessenbergQ{<:Any,false}) where {T} = convert(Matrix{T}, LAPACK.orghr!(1, size(Q.factors, 1), copy(Q.factors), Q.τ)) -Matrix{T}(Q::BlasHessenbergQ{<:Any,true}) where {T} = convert(Matrix{T}, LAPACK.orgtr!(Q.uplo, copy(Q.factors), Q.τ)) AbstractArray(F::Hessenberg) = AbstractMatrix(F) Matrix(F::Hessenberg) = Array(AbstractArray(F)) Array(F::Hessenberg) = Matrix(F) @@ -574,31 +519,6 @@ 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::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::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::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::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,<: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,<: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 function rmul_triu!(M::AbstractMatrix, x, k::Integer=0) diff --git a/stdlib/LinearAlgebra/src/lq.jl b/stdlib/LinearAlgebra/src/lq.jl index a162333be339a..33d794906c7e6 100644 --- a/stdlib/LinearAlgebra/src/lq.jl +++ b/stdlib/LinearAlgebra/src/lq.jl @@ -28,9 +28,7 @@ L factor: -8.60233 0.0 4.41741 -0.697486 Q factor: -2×2 LinearAlgebra.LQPackedQ{Float64, Matrix{Float64}, Vector{Float64}}: - -0.581238 -0.813733 - -0.813733 0.581238 +2×2 LinearAlgebra.LQPackedQ{Float64, Matrix{Float64}, Vector{Float64}} julia> S.L * S.Q 2×2 Matrix{Float64}: @@ -65,12 +63,6 @@ Base.iterate(S::LQ) = (S.L, Val(:Q)) Base.iterate(S::LQ, ::Val{:Q}) = (S.Q, Val(:done)) Base.iterate(S::LQ, ::Val{:done}) = nothing -struct LQPackedQ{T,S<:AbstractMatrix{T},C<:AbstractVector{T}} <: AbstractMatrix{T} - factors::S - τ::C -end - - """ lq!(A) -> LQ @@ -78,6 +70,7 @@ Compute the [`LQ`](@ref) factorization of `A`, using the input matrix as a workspace. See also [`lq`](@ref). """ lq!(A::StridedMatrix{<:BlasFloat}) = LQ(LAPACK.gelqf!(A)...) + """ lq(A) -> S::LQ @@ -105,9 +98,7 @@ L factor: -8.60233 0.0 4.41741 -0.697486 Q factor: -2×2 LinearAlgebra.LQPackedQ{Float64, Matrix{Float64}, Vector{Float64}}: - -0.581238 -0.813733 - -0.813733 0.581238 +2×2 LinearAlgebra.LQPackedQ{Float64, Matrix{Float64}, Vector{Float64}} julia> S.L * S.Q 2×2 Matrix{Float64}: @@ -156,8 +147,8 @@ end Base.propertynames(F::LQ, private::Bool=false) = (:L, :Q, (private ? fieldnames(typeof(F)) : ())...) -getindex(A::LQPackedQ, i::Integer, j::Integer) = - lmul!(A, setindex!(zeros(eltype(A), size(A, 2)), 1, j))[i] +# getindex(A::LQPackedQ, i::Integer, j::Integer) = +# lmul!(A, setindex!(zeros(eltype(A), size(A, 2)), 1, j))[i] function show(io::IO, mime::MIME{Symbol("text/plain")}, F::LQ) summary(io, F); println(io) @@ -167,32 +158,9 @@ function show(io::IO, mime::MIME{Symbol("text/plain")}, F::LQ) show(io, mime, F.Q) end -LQPackedQ{T}(Q::LQPackedQ) where {T} = LQPackedQ(convert(AbstractMatrix{T}, Q.factors), convert(Vector{T}, Q.τ)) -AbstractMatrix{T}(Q::LQPackedQ) where {T} = LQPackedQ{T}(Q) -Matrix{T}(A::LQPackedQ) where {T} = convert(Matrix{T}, LAPACK.orglq!(copy(A.factors),A.τ)) -Matrix(A::LQPackedQ{T}) where {T} = Matrix{T}(A) -Array{T}(A::LQPackedQ{T}) where {T} = Matrix{T}(A) -Array(A::LQPackedQ) = Matrix(A) - size(F::LQ, dim::Integer) = size(getfield(F, :factors), dim) size(F::LQ) = size(getfield(F, :factors)) -# size(Q::LQPackedQ) yields the shape of Q's square form -function size(Q::LQPackedQ) - n = size(Q.factors, 2) - return n, n -end -function size(Q::LQPackedQ, dim::Integer) - if dim < 1 - throw(BoundsError()) - elseif dim <= 2 # && 1 <= dim - return size(Q.factors, 2) - else # 2 < dim - return 1 - end -end - - ## Multiplication by LQ function lmul!(A::LQ, B::AbstractVecOrMat) lmul!(LowerTriangular(A.L), view(lmul!(A.Q, B), 1:size(A,1), axes(B,2))) @@ -203,127 +171,6 @@ function *(A::LQ{TA}, B::AbstractVecOrMat{TB}) where {TA,TB} _cut_B(lmul!(convert(Factorization{TAB}, A), copy_similar(B, TAB)), 1:size(A,1)) end -## Multiplication by Q -### QB -lmul!(A::LQPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = LAPACK.ormlq!('L','N',A.factors,A.τ,B) -function (*)(A::LQPackedQ, B::StridedVecOrMat) - TAB = promote_type(eltype(A), eltype(B)) - lmul!(AbstractMatrix{TAB}(A), copymutable_oftype(B, TAB)) -end - -### QcB -lmul!(adjA::Adjoint{<:Any,<:LQPackedQ{T}}, B::StridedVecOrMat{T}) where {T<:BlasReal} = - (A = adjA.parent; LAPACK.ormlq!('L', 'T', A.factors, A.τ, B)) -lmul!(adjA::Adjoint{<:Any,<:LQPackedQ{T}}, B::StridedVecOrMat{T}) where {T<:BlasComplex} = - (A = adjA.parent; LAPACK.ormlq!('L', 'C', A.factors, A.τ, B)) - -function *(adjA::Adjoint{<:Any,<:LQPackedQ}, B::StridedVecOrMat) - A = adjA.parent - TAB = promote_type(eltype(A), eltype(B)) - if size(B,1) == size(A.factors,2) - lmul!(adjoint(AbstractMatrix{TAB}(A)), copymutable_oftype(B, TAB)) - elseif size(B,1) == size(A.factors,1) - lmul!(adjoint(AbstractMatrix{TAB}(A)), [B; zeros(TAB, size(A.factors, 2) - size(A.factors, 1), size(B, 2))]) - else - throw(DimensionMismatch("first dimension of B, $(size(B,1)), must equal one of the dimensions of A, $(size(A))")) - end -end - -### QBc/QcBc -function *(A::LQPackedQ, adjB::Adjoint{<:Any,<:StridedVecOrMat}) - B = adjB.parent - TAB = promote_type(eltype(A), eltype(B)) - BB = similar(B, TAB, (size(B, 2), size(B, 1))) - adjoint!(BB, B) - return lmul!(A, BB) -end -function *(adjA::Adjoint{<:Any,<:LQPackedQ}, adjB::Adjoint{<:Any,<:StridedVecOrMat}) - B = adjB.parent - TAB = promote_type(eltype(adjA.parent), eltype(B)) - BB = similar(B, TAB, (size(B, 2), size(B, 1))) - adjoint!(BB, B) - return lmul!(adjA, BB) -end - -# in-place right-application of LQPackedQs -# these methods require that the applied-to matrix's (A's) number of columns -# match the number of columns (nQ) of the LQPackedQ (Q) (necessary for in-place -# operation, and the underlying LAPACK routine (ormlq) treats the implicit Q -# as its (nQ-by-nQ) square form) -rmul!(A::StridedMatrix{T}, B::LQPackedQ{T}) where {T<:BlasFloat} = - LAPACK.ormlq!('R', 'N', B.factors, B.τ, A) -rmul!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:LQPackedQ{T}}) where {T<:BlasReal} = - (B = adjB.parent; LAPACK.ormlq!('R', 'T', B.factors, B.τ, A)) -rmul!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:LQPackedQ{T}}) where {T<:BlasComplex} = - (B = adjB.parent; LAPACK.ormlq!('R', 'C', B.factors, B.τ, A)) - -# out-of-place right application of LQPackedQs -# -# LQPackedQ's out-of-place multiplication behavior is context dependent. specifically, -# if the inner dimension in the multiplication is the LQPackedQ's second dimension, -# the LQPackedQ behaves like its square form. if the inner dimension in the -# multiplication is the LQPackedQ's first dimension, the LQPackedQ behaves like either -# its square form or its truncated form depending on the shape of the other object -# involved in the multiplication. we treat these cases separately. -# -# (1) the inner dimension in the multiplication is the LQPackedQ's second dimension. -# in this case, the LQPackedQ behaves like its square form. -# -function *(A::StridedVecOrMat, adjQ::Adjoint{<:Any,<:LQPackedQ}) - Q = adjQ.parent - TR = promote_type(eltype(A), eltype(Q)) - return rmul!(copymutable_oftype(A, TR), adjoint(AbstractMatrix{TR}(Q))) -end -function *(adjA::Adjoint{<:Any,<:StridedMatrix}, adjQ::Adjoint{<:Any,<:LQPackedQ}) - A, Q = adjA.parent, adjQ.parent - TR = promote_type(eltype(A), eltype(Q)) - C = adjoint!(similar(A, TR, reverse(size(A))), A) - return rmul!(C, adjoint(AbstractMatrix{TR}(Q))) -end -# -# (2) the inner dimension in the multiplication is the LQPackedQ's first dimension. -# in this case, the LQPackedQ behaves like either its square form or its -# truncated form depending on the shape of the other object in the multiplication. -# -# these methods: (1) check whether the applied-to matrix's (A's) appropriate dimension -# (columns for A_*, rows for Ac_*) matches the number of columns (nQ) of the LQPackedQ (Q), -# and if so effectively apply Q's square form to A without additional shenanigans; and -# (2) if the preceding dimensions do not match, check whether the appropriate dimension of -# A instead matches the number of rows of the matrix of which Q is a factor (i.e. -# size(Q.factors, 1)), and if so implicitly apply Q's truncated form to A by zero extending -# A as necessary for check (1) to pass (if possible) and then applying Q's square form -# -function *(A::StridedVecOrMat, Q::LQPackedQ) - TR = promote_type(eltype(A), eltype(Q)) - if size(A, 2) == size(Q.factors, 2) - C = copymutable_oftype(A, TR) - elseif size(A, 2) == size(Q.factors, 1) - C = zeros(TR, size(A, 1), size(Q.factors, 2)) - copyto!(C, 1, A, 1, length(A)) - else - _rightappdimmismatch("columns") - end - return rmul!(C, AbstractMatrix{TR}(Q)) -end -function *(adjA::Adjoint{<:Any,<:StridedMatrix}, Q::LQPackedQ) - A = adjA.parent - TR = promote_type(eltype(A), eltype(Q)) - if size(A, 1) == size(Q.factors, 2) - C = adjoint!(similar(A, TR, reverse(size(A))), A) - elseif size(A, 1) == size(Q.factors, 1) - C = zeros(TR, size(A, 2), size(Q.factors, 2)) - adjoint!(view(C, :, 1:size(A, 1)), A) - else - _rightappdimmismatch("rows") - end - return rmul!(C, AbstractMatrix{TR}(Q)) -end -_rightappdimmismatch(rowsorcols) = - throw(DimensionMismatch(string("the number of $(rowsorcols) of the matrix on the left ", - "must match either (1) the number of columns of the (LQPackedQ) matrix on the right ", - "or (2) the number of rows of that (LQPackedQ) matrix's internal representation ", - "(the factorization's originating matrix's number of rows)"))) - # 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 @@ -356,7 +203,3 @@ function ldiv!(Fadj::AdjointFactorization{<:Any,<:LQ}, B::AbstractVecOrMat) 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. -det(Q::LQPackedQ) = conj(_det_tau(Q.τ)) diff --git a/stdlib/LinearAlgebra/src/qr.jl b/stdlib/LinearAlgebra/src/qr.jl index 30deb785e6019..d2420fb78edef 100644 --- a/stdlib/LinearAlgebra/src/qr.jl +++ b/stdlib/LinearAlgebra/src/qr.jl @@ -32,7 +32,6 @@ The object has two fields: ``v_i`` is the ``i``th column of the matrix `V = I + tril(F.factors, -1)`. * `τ` is a vector of length `min(m,n)` containing the coefficients ``\tau_i``. - """ struct QR{T,S<:AbstractMatrix{T},C<:AbstractVector{T}} <: Factorization{T} factors::S @@ -298,7 +297,7 @@ qr!(A::StridedMatrix{<:BlasFloat}, ::ColumnNorm) = QRPivoted(LAPACK.geqp3!(A)... """ qr!(A, pivot = NoPivot(); blocksize) -`qr!` is the same as [`qr`](@ref) when `A` is a subtype of [`StridedMatrix`](@ref), +`qr!` is the same as [`qr`](@ref) when `A` is a subtype of [`AbstractMatrix`](@ref), but saves space by overwriting the input `A`, instead of creating a copy. An [`InexactError`](@ref) exception is thrown if the factorization produces a number not representable by the element type of `A`, e.g. for integer types. @@ -316,9 +315,7 @@ julia> a = [1. 2.; 3. 4.] julia> qr!(a) LinearAlgebra.QRCompactWY{Float64, Matrix{Float64}, Matrix{Float64}} Q factor: -2×2 LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}, Matrix{Float64}}: - -0.316228 -0.948683 - -0.948683 0.316228 +2×2 LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}, Matrix{Float64}} R factor: 2×2 Matrix{Float64}: -3.16228 -4.42719 @@ -387,7 +384,7 @@ orthogonal matrix. The block size for QR decomposition can be specified by keyword argument `blocksize :: Integer` when `pivot == NoPivot()` and `A isa StridedMatrix{<:BlasFloat}`. -It is ignored when `blocksize > minimum(size(A))`. See [`QRCompactWY`](@ref). +It is ignored when `blocksize > minimum(size(A))`. See [`QRCompactWY`](@ref). !!! compat "Julia 1.4" The `blocksize` keyword argument requires Julia 1.4 or later. @@ -403,10 +400,7 @@ julia> A = [3.0 -6.0; 4.0 -8.0; 0.0 1.0] julia> F = qr(A) LinearAlgebra.QRCompactWY{Float64, Matrix{Float64}, Matrix{Float64}} Q factor: -3×3 LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}, Matrix{Float64}}: - -0.6 0.0 0.8 - -0.8 0.0 -0.6 - 0.0 -1.0 0.0 +3×3 LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}, Matrix{Float64}} R factor: 2×2 Matrix{Float64}: -5.0 10.0 @@ -518,372 +512,25 @@ transpose(F::Union{QR{<:Real},QRPivoted{<:Real},QRCompactWY{<:Real}}) = F' transpose(::Union{QR,QRPivoted,QRCompactWY}) = throw(ArgumentError("transpose of QR decomposition is not supported, consider using adjoint")) -abstract type AbstractQ{T} <: AbstractMatrix{T} end - -inv(Q::AbstractQ) = Q' - -""" - QRPackedQ <: AbstractMatrix - -The orthogonal/unitary ``Q`` matrix of a QR factorization stored in [`QR`](@ref) or -[`QRPivoted`](@ref) format. -""" -struct QRPackedQ{T,S<:AbstractMatrix{T},C<:AbstractVector{T}} <: AbstractQ{T} - factors::S - τ::C - - function QRPackedQ{T,S,C}(factors, τ) where {T,S<:AbstractMatrix{T},C<:AbstractVector{T}} - require_one_based_indexing(factors) - new{T,S,C}(factors, τ) - end -end -QRPackedQ(factors::AbstractMatrix{T}, τ::AbstractVector{T}) where {T} = - QRPackedQ{T,typeof(factors),typeof(τ)}(factors, τ) -QRPackedQ{T}(factors::AbstractMatrix, τ::AbstractVector) where {T} = - QRPackedQ(convert(AbstractMatrix{T}, factors), convert(AbstractVector{T}, τ)) -# backwards-compatible constructors (remove with Julia 2.0) -@deprecate(QRPackedQ{T,S}(factors::AbstractMatrix{T}, τ::AbstractVector{T}) where {T,S}, - QRPackedQ{T,S,typeof(τ)}(factors, τ), false) - -""" - QRCompactWYQ <: AbstractMatrix - -The orthogonal/unitary ``Q`` matrix of a QR factorization stored in [`QRCompactWY`](@ref) -format. -""" -struct QRCompactWYQ{S, M<:AbstractMatrix{S}, C<:AbstractMatrix{S}} <: AbstractQ{S} - factors::M - T::C - - function QRCompactWYQ{S,M,C}(factors, T) where {S,M<:AbstractMatrix{S},C<:AbstractMatrix{S}} - require_one_based_indexing(factors) - new{S,M,C}(factors, T) - end -end -QRCompactWYQ(factors::AbstractMatrix{S}, T::AbstractMatrix{S}) where {S} = - QRCompactWYQ{S,typeof(factors),typeof(T)}(factors, T) -QRCompactWYQ{S}(factors::AbstractMatrix, T::AbstractMatrix) where {S} = - QRCompactWYQ(convert(AbstractMatrix{S}, factors), convert(AbstractMatrix{S}, T)) -# backwards-compatible constructors (remove with Julia 2.0) -@deprecate(QRCompactWYQ{S,M}(factors::AbstractMatrix{S}, T::AbstractMatrix{S}) where {S,M}, - QRCompactWYQ{S,M,typeof(T)}(factors, T), false) - -QRPackedQ{T}(Q::QRPackedQ) where {T} = QRPackedQ(convert(AbstractMatrix{T}, Q.factors), convert(Vector{T}, Q.τ)) -AbstractMatrix{T}(Q::QRPackedQ{T}) where {T} = Q -AbstractMatrix{T}(Q::QRPackedQ) where {T} = QRPackedQ{T}(Q) -QRCompactWYQ{S}(Q::QRCompactWYQ) where {S} = QRCompactWYQ(convert(AbstractMatrix{S}, Q.factors), convert(AbstractMatrix{S}, Q.T)) -AbstractMatrix{S}(Q::QRCompactWYQ{S}) where {S} = Q -AbstractMatrix{S}(Q::QRCompactWYQ) where {S} = QRCompactWYQ{S}(Q) -Matrix{T}(Q::AbstractQ{S}) where {T,S} = convert(Matrix{T}, lmul!(Q, Matrix{S}(I, size(Q, 1), min(size(Q.factors)...)))) -Matrix(Q::AbstractQ{T}) where {T} = Matrix{T}(Q) -Array{T}(Q::AbstractQ) where {T} = Matrix{T}(Q) -Array(Q::AbstractQ) = Matrix(Q) - -size(F::Union{QR,QRCompactWY,QRPivoted}, dim::Integer) = size(getfield(F, :factors), dim) size(F::Union{QR,QRCompactWY,QRPivoted}) = size(getfield(F, :factors)) -size(Q::Union{QRCompactWYQ,QRPackedQ}, dim::Integer) = - size(getfield(Q, :factors), dim == 2 ? 1 : dim) -size(Q::Union{QRCompactWYQ,QRPackedQ}) = size(Q, 1), size(Q, 2) - -copymutable(Q::AbstractQ{T}) where {T} = lmul!(Q, Matrix{T}(I, size(Q))) -copy(Q::AbstractQ) = copymutable(Q) -getindex(Q::AbstractQ, inds...) = copymutable(Q)[inds...] -getindex(Q::AbstractQ, ::Colon, ::Colon) = copy(Q) - -function getindex(Q::AbstractQ, ::Colon, j::Int) - y = zeros(eltype(Q), size(Q, 2)) - y[j] = 1 - lmul!(Q, y) -end - -getindex(Q::AbstractQ, i::Int, j::Int) = Q[:, j][i] - -# specialization avoiding the fallback using slow `getindex` -function copyto!(dest::AbstractMatrix, src::AbstractQ) - copyto!(dest, I) - lmul!(src, dest) -end -# needed to resolve method ambiguities -function copyto!(dest::PermutedDimsArray{T,2,perm}, src::AbstractQ) where {T,perm} - if perm == (1, 2) - copyto!(parent(dest), src) - else - @assert perm == (2, 1) # there are no other permutations of two indices - if T <: Real - copyto!(parent(dest), I) - lmul!(src', parent(dest)) - else - # LAPACK does not offer inplace lmul!(transpose(Q), B) for complex Q - tmp = similar(parent(dest)) - copyto!(tmp, I) - rmul!(tmp, src) - permutedims!(parent(dest), tmp, (2, 1)) - end - end - return dest -end - -## Multiplication by Q -### QB -lmul!(A::QRCompactWYQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasFloat, S<:StridedMatrix} = - LAPACK.gemqrt!('L', 'N', A.factors, A.T, B) -lmul!(A::QRPackedQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasFloat, S<:StridedMatrix} = - LAPACK.ormqr!('L', 'N', A.factors, A.τ, B) -function lmul!(A::QRPackedQ, B::AbstractVecOrMat) - require_one_based_indexing(B) - mA, nA = size(A.factors) - mB, nB = size(B,1), size(B,2) - if mA != mB - throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)")) - end - Afactors = A.factors - @inbounds begin - for k = min(mA,nA):-1:1 - for j = 1:nB - vBj = B[k,j] - for i = k+1:mB - vBj += conj(Afactors[i,k])*B[i,j] - end - vBj = A.τ[k]*vBj - B[k,j] -= vBj - for i = k+1:mB - B[i,j] -= Afactors[i,k]*vBj - end - end - end - end - B -end - -function (*)(A::AbstractQ, b::StridedVector) - TAb = promote_type(eltype(A), eltype(b)) - Anew = convert(AbstractMatrix{TAb}, A) - if size(A.factors, 1) == length(b) - bnew = copymutable_oftype(b, TAb) - elseif size(A.factors, 2) == length(b) - bnew = [b; zeros(TAb, size(A.factors, 1) - length(b))] - else - throw(DimensionMismatch("vector must have length either $(size(A.factors, 1)) or $(size(A.factors, 2))")) - end - lmul!(Anew, bnew) -end -function (*)(A::AbstractQ, B::StridedMatrix) - TAB = promote_type(eltype(A), eltype(B)) - Anew = convert(AbstractMatrix{TAB}, A) - if size(A.factors, 1) == size(B, 1) - Bnew = copymutable_oftype(B, TAB) - elseif size(A.factors, 2) == size(B, 1) - Bnew = [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 - -function (*)(A::AbstractQ, b::Number) - TAb = promote_type(eltype(A), typeof(b)) - dest = similar(A, TAb) - copyto!(dest, b*I) - lmul!(A, dest) -end - -### QcB -lmul!(adjA::Adjoint{<:Any,<:QRCompactWYQ{T,S}}, B::StridedVecOrMat{T}) where {T<:BlasReal,S<:StridedMatrix} = - (A = adjA.parent; LAPACK.gemqrt!('L', 'T', A.factors, A.T, B)) -lmul!(adjA::Adjoint{<:Any,<:QRCompactWYQ{T,S}}, B::StridedVecOrMat{T}) where {T<:BlasComplex,S<:StridedMatrix} = - (A = adjA.parent; LAPACK.gemqrt!('L', 'C', A.factors, A.T, B)) -lmul!(adjA::Adjoint{<:Any,<:QRPackedQ{T,S}}, B::StridedVecOrMat{T}) where {T<:BlasReal,S<:StridedMatrix} = - (A = adjA.parent; LAPACK.ormqr!('L', 'T', A.factors, A.τ, B)) -lmul!(adjA::Adjoint{<:Any,<:QRPackedQ{T,S}}, B::StridedVecOrMat{T}) where {T<:BlasComplex,S<:StridedMatrix} = - (A = adjA.parent; LAPACK.ormqr!('L', 'C', A.factors, A.τ, B)) -function lmul!(adjA::Adjoint{<:Any,<:QRPackedQ}, B::AbstractVecOrMat) - require_one_based_indexing(B) - A = adjA.parent - mA, nA = size(A.factors) - mB, nB = size(B,1), size(B,2) - if mA != mB - throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)")) - end - Afactors = A.factors - @inbounds begin - for k = 1:min(mA,nA) - for j = 1:nB - vBj = B[k,j] - for i = k+1:mB - vBj += conj(Afactors[i,k])*B[i,j] - end - vBj = conj(A.τ[k])*vBj - B[k,j] -= vBj - for i = k+1:mB - B[i,j] -= Afactors[i,k]*vBj - end - end - end - end - B -end -function *(adjQ::Adjoint{<:Any,<:AbstractQ}, B::StridedVecOrMat) - Q = adjQ.parent - TQB = promote_type(eltype(Q), eltype(B)) - return lmul!(adjoint(convert(AbstractMatrix{TQB}, Q)), copymutable_oftype(B, TQB)) -end - -### QBc/QcBc -function *(Q::AbstractQ, adjB::Adjoint{<:Any,<:StridedVecOrMat}) - B = adjB.parent - TQB = promote_type(eltype(Q), eltype(B)) - Bc = similar(B, TQB, (size(B, 2), size(B, 1))) - adjoint!(Bc, B) - return lmul!(convert(AbstractMatrix{TQB}, Q), Bc) -end -function *(adjQ::Adjoint{<:Any,<:AbstractQ}, adjB::Adjoint{<:Any,<:StridedVecOrMat}) - Q, B = adjQ.parent, adjB.parent - TQB = promote_type(eltype(Q), eltype(B)) - Bc = similar(B, TQB, (size(B, 2), size(B, 1))) - adjoint!(Bc, B) - return lmul!(adjoint(convert(AbstractMatrix{TQB}, Q)), Bc) -end - -### AQ -rmul!(A::StridedVecOrMat{T}, B::QRCompactWYQ{T,S}) where {T<:BlasFloat,S<:StridedMatrix} = - LAPACK.gemqrt!('R', 'N', B.factors, B.T, A) -rmul!(A::StridedVecOrMat{T}, B::QRPackedQ{T,S}) where {T<:BlasFloat,S<:StridedMatrix} = - LAPACK.ormqr!('R', 'N', B.factors, B.τ, A) -function rmul!(A::StridedMatrix,Q::QRPackedQ) - mQ, nQ = size(Q.factors) - mA, nA = size(A,1), size(A,2) - if nA != mQ - throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but matrix Q has dimensions ($mQ, $nQ)")) - end - Qfactors = Q.factors - @inbounds begin - for k = 1:min(mQ,nQ) - for i = 1:mA - vAi = A[i,k] - for j = k+1:mQ - vAi += A[i,j]*Qfactors[j,k] - end - vAi = vAi*Q.τ[k] - A[i,k] -= vAi - for j = k+1:nA - A[i,j] -= vAi*conj(Qfactors[j,k]) - end - end - end - end - A -end - -function (*)(A::StridedMatrix, Q::AbstractQ) - TAQ = promote_type(eltype(A), eltype(Q)) - - return rmul!(copymutable_oftype(A, TAQ), convert(AbstractMatrix{TAQ}, Q)) -end - -function (*)(a::Number, B::AbstractQ) - TaB = promote_type(typeof(a), eltype(B)) - dest = similar(B, TaB) - copyto!(dest, a*I) - rmul!(dest, B) -end - -### AQc -rmul!(A::StridedVecOrMat{T}, adjB::Adjoint{<:Any,<:QRCompactWYQ{T}}) where {T<:BlasReal} = - (B = adjB.parent; LAPACK.gemqrt!('R', 'T', B.factors, B.T, A)) -rmul!(A::StridedVecOrMat{T}, adjB::Adjoint{<:Any,<:QRCompactWYQ{T}}) where {T<:BlasComplex} = - (B = adjB.parent; LAPACK.gemqrt!('R', 'C', B.factors, B.T, A)) -rmul!(A::StridedVecOrMat{T}, adjB::Adjoint{<:Any,<:QRPackedQ{T}}) where {T<:BlasReal} = - (B = adjB.parent; LAPACK.ormqr!('R', 'T', B.factors, B.τ, A)) -rmul!(A::StridedVecOrMat{T}, adjB::Adjoint{<:Any,<:QRPackedQ{T}}) where {T<:BlasComplex} = - (B = adjB.parent; LAPACK.ormqr!('R', 'C', B.factors, B.τ, A)) -function rmul!(A::StridedMatrix, adjQ::Adjoint{<:Any,<:QRPackedQ}) - Q = adjQ.parent - mQ, nQ = size(Q.factors) - mA, nA = size(A,1), size(A,2) - if nA != mQ - throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but matrix Q has dimensions ($mQ, $nQ)")) - end - Qfactors = Q.factors - @inbounds begin - for k = min(mQ,nQ):-1:1 - for i = 1:mA - vAi = A[i,k] - for j = k+1:mQ - vAi += A[i,j]*Qfactors[j,k] - end - vAi = vAi*conj(Q.τ[k]) - A[i,k] -= vAi - for j = k+1:nA - A[i,j] -= vAi*conj(Qfactors[j,k]) - end - end - end - end - A -end -function *(A::StridedMatrix, adjB::Adjoint{<:Any,<:AbstractQ}) - B = adjB.parent - TAB = promote_type(eltype(A),eltype(B)) - BB = convert(AbstractMatrix{TAB}, B) - if size(A,2) == size(B.factors, 1) - AA = copy_similar(A, TAB) - return rmul!(AA, adjoint(BB)) - elseif size(A,2) == size(B.factors,2) - return rmul!([A zeros(TAB, size(A, 1), size(B.factors, 1) - size(B.factors, 2))], adjoint(BB)) - else - throw(DimensionMismatch("matrix A has dimensions $(size(A)) but matrix B has dimensions $(size(B))")) - end -end -*(u::AdjointAbsVec, A::Adjoint{<:Any,<:AbstractQ}) = adjoint(A.parent * u.parent) - - -### AcQ/AcQc -function *(adjA::Adjoint{<:Any,<:StridedVecOrMat}, Q::AbstractQ) - A = adjA.parent - TAQ = promote_type(eltype(A), eltype(Q)) - Ac = similar(A, TAQ, (size(A, 2), size(A, 1))) - adjoint!(Ac, A) - return rmul!(Ac, convert(AbstractMatrix{TAQ}, Q)) -end -function *(adjA::Adjoint{<:Any,<:StridedVecOrMat}, adjQ::Adjoint{<:Any,<:AbstractQ}) - A, Q = adjA.parent, adjQ.parent - TAQ = promote_type(eltype(A), eltype(Q)) - Ac = similar(A, TAQ, (size(A, 2), size(A, 1))) - adjoint!(Ac, A) - return rmul!(Ac, adjoint(convert(AbstractMatrix{TAQ}, Q))) -end +size(F::Union{QR,QRCompactWY,QRPivoted}, dim::Integer) = size(getfield(F, :factors), dim) -### mul! -function mul!(C::StridedVecOrMat{T}, Q::AbstractQ{T}, B::StridedVecOrMat{T}) where {T} - require_one_based_indexing(C, B) - mB = size(B, 1) - mC = size(C, 1) - if mB < mC - inds = CartesianIndices(B) - copyto!(C, inds, B, inds) - C[CartesianIndices((mB+1:mC, axes(C, 2)))] .= zero(T) - return lmul!(Q, C) - else - return lmul!(Q, copyto!(C, B)) - end -end -mul!(C::StridedVecOrMat{T}, A::StridedVecOrMat{T}, Q::AbstractQ{T}) where {T} = rmul!(copyto!(C, A), Q) -mul!(C::StridedVecOrMat{T}, adjQ::Adjoint{<:Any,<:AbstractQ{T}}, B::StridedVecOrMat{T}) where {T} = lmul!(adjQ, copyto!(C, B)) -mul!(C::StridedVecOrMat{T}, A::StridedVecOrMat{T}, adjQ::Adjoint{<:Any,<:AbstractQ{T}}) where {T} = rmul!(copyto!(C, A), adjQ) -function ldiv!(A::QRCompactWY{T}, b::AbstractVector{T}) where {T<:BlasFloat} - m,n = size(A) +function ldiv!(A::QRCompactWY{T}, b::AbstractVector{T}) where {T} + require_one_based_indexing(b) + m, n = size(A) ldiv!(UpperTriangular(view(A.factors, 1:min(m,n), 1:n)), view(lmul!(adjoint(A.Q), b), 1:size(A, 2))) return b end -function ldiv!(A::QRCompactWY{T}, B::AbstractMatrix{T}) where {T<:BlasFloat} - m,n = size(A) +function ldiv!(A::QRCompactWY{T}, B::AbstractMatrix{T}) where {T} + require_one_based_indexing(B) + m, n = size(A) ldiv!(UpperTriangular(view(A.factors, 1:min(m,n), 1:n)), view(lmul!(adjoint(A.Q), B), 1:size(A, 2), 1:size(B, 2))) return B end # Julia implementation similar to xgelsy -function ldiv!(A::QRPivoted{T}, B::AbstractMatrix{T}, rcond::Real) where T<:BlasFloat +function ldiv!(A::QRPivoted{T,<:StridedMatrix}, B::AbstractMatrix{T}, rcond::Real) where {T<:BlasFloat} mA, nA = size(A.factors) nr = min(mA,nA) nrhs = size(B, 2) @@ -919,9 +566,9 @@ function ldiv!(A::QRPivoted{T}, B::AbstractMatrix{T}, rcond::Real) where T<:Blas B[A.p,:] = B[1:nA,:] return B, rnk end -ldiv!(A::QRPivoted{T}, B::AbstractVector{T}) where {T<:BlasFloat} = +ldiv!(A::QRPivoted{T,<:StridedMatrix}, B::AbstractVector{T}) where {T<:BlasFloat} = vec(ldiv!(A, reshape(B, length(B), 1))) -ldiv!(A::QRPivoted{T}, B::AbstractVecOrMat{T}) where {T<:BlasFloat} = +ldiv!(A::QRPivoted{T,<:StridedMatrix}, B::AbstractMatrix{T}) where {T<:BlasFloat} = ldiv!(A, B, min(size(A)...)*eps(real(float(one(eltype(B))))))[1] function _wide_qr_ldiv!(A::QR{T}, B::AbstractMatrix{T}) where T m, n = size(A) @@ -1002,7 +649,7 @@ function _apply_permutation!(F::QRPivoted, B::AbstractVecOrMat) B[1:length(F.p), :] = B[F.p, :] return B end -_apply_permutation!(F::Factorization, B::AbstractVecOrMat) = B +_apply_permutation!(::Factorization, B::AbstractVecOrMat) = B function ldiv!(Fadj::AdjointFactorization{<:Any,<:Union{QR,QRCompactWY,QRPivoted}}, B::AbstractVecOrMat) require_one_based_indexing(B) @@ -1063,25 +710,3 @@ end ## Lower priority: Add LQ, QL and RQ factorizations # FIXME! Should add balancing option through xgebal - - -det(Q::QRPackedQ) = _det_tau(Q.τ) - -det(Q::QRCompactWYQ) = - prod(i -> _det_tau(_diagview(Q.T[:, i:min(i + size(Q.T, 1), size(Q.T, 2))])), - 1:size(Q.T, 1):size(Q.T, 2)) - -_diagview(A) = @view A[diagind(A)] - -# Compute `det` from the number of Householder reflections. Handle -# the case `Q.τ` contains zeros. -_det_tau(τs::AbstractVector{<:Real}) = - isodd(count(!iszero, τs)) ? -one(eltype(τs)) : one(eltype(τs)) - -# In complex case, we need to compute the non-unit eigenvalue `λ = 1 - c*τ` -# (where `c = v'v`) of each Householder reflector. As we know that the -# reflector must have the determinant of 1, it must satisfy `abs2(λ) == 1`. -# Combining this with the constraint `c > 0`, it turns out that the eigenvalue -# (hence the determinant) can be computed as `λ = -sign(τ)^2`. -# See: https://github.com/JuliaLang/julia/pull/32887#issuecomment-521935716 -_det_tau(τs) = prod(τ -> iszero(τ) ? one(τ) : -sign(τ)^2, τs) diff --git a/stdlib/LinearAlgebra/src/special.jl b/stdlib/LinearAlgebra/src/special.jl index e4f28286b6aaa..19b1057f9d6b8 100644 --- a/stdlib/LinearAlgebra/src/special.jl +++ b/stdlib/LinearAlgebra/src/special.jl @@ -107,6 +107,29 @@ for op in (:+, :-) end end +function *(H::UpperHessenberg, B::Bidiagonal) + T = promote_op(matprod, eltype(H), eltype(B)) + A = mul!(similar(H, T, size(H)), H, B) + return B.uplo == 'U' ? UpperHessenberg(A) : A +end +function *(B::Bidiagonal, H::UpperHessenberg) + T = promote_op(matprod, eltype(B), eltype(H)) + A = mul!(similar(H, T, size(H)), B, H) + return B.uplo == 'U' ? UpperHessenberg(A) : A +end + +function /(H::UpperHessenberg, B::Bidiagonal) + T = typeof(oneunit(eltype(H))/oneunit(eltype(B))) + A = _rdiv!(similar(H, T, size(H)), H, B) + return B.uplo == 'U' ? UpperHessenberg(A) : A +end + +function \(B::Bidiagonal, H::UpperHessenberg) + T = typeof(oneunit(eltype(B))\oneunit(eltype(H))) + A = ldiv!(similar(H, T, size(H)), B, H) + return B.uplo == 'U' ? UpperHessenberg(A) : A +end + # specialized +/- for structured matrices. If these are removed, it falls # back to broadcasting which has ~2-10x speed regressions. # For the other structure matrix pairs, broadcasting works well. @@ -235,65 +258,15 @@ function (-)(A::UniformScaling, B::Diagonal) Diagonal(Ref(A) .- B.diag) end -lmul!(Q::AbstractQ, B::AbstractTriangular) = lmul!(Q, full!(B)) -lmul!(Q::QRPackedQ, B::AbstractTriangular) = lmul!(Q, full!(B)) # disambiguation -lmul!(Q::Adjoint{<:Any,<:AbstractQ}, B::AbstractTriangular) = lmul!(Q, full!(B)) -lmul!(Q::Adjoint{<:Any,<:QRPackedQ}, B::AbstractTriangular) = lmul!(Q, full!(B)) # disambiguation - -function _qlmul(Q::AbstractQ, B) - TQB = promote_type(eltype(Q), eltype(B)) - if size(Q.factors, 1) == size(B, 1) - Bnew = Matrix{TQB}(B) - elseif size(Q.factors, 2) == size(B, 1) - Bnew = [Matrix{TQB}(B); zeros(TQB, size(Q.factors, 1) - size(B,1), size(B, 2))] - else - throw(DimensionMismatch("first dimension of matrix must have size either $(size(Q.factors, 1)) or $(size(Q.factors, 2))")) - end - lmul!(convert(AbstractMatrix{TQB}, Q), Bnew) -end -function _qlmul(adjQ::Adjoint{<:Any,<:AbstractQ}, B) - TQB = promote_type(eltype(adjQ), eltype(B)) - lmul!(adjoint(convert(AbstractMatrix{TQB}, parent(adjQ))), Matrix{TQB}(B)) -end - -*(Q::AbstractQ, B::AbstractTriangular) = _qlmul(Q, B) -*(Q::Adjoint{<:Any,<:AbstractQ}, B::AbstractTriangular) = _qlmul(Q, B) -*(Q::AbstractQ, B::BiTriSym) = _qlmul(Q, B) -*(Q::Adjoint{<:Any,<:AbstractQ}, B::BiTriSym) = _qlmul(Q, B) -*(Q::AbstractQ, B::Diagonal) = _qlmul(Q, B) -*(Q::Adjoint{<:Any,<:AbstractQ}, B::Diagonal) = _qlmul(Q, B) - -rmul!(A::AbstractTriangular, Q::AbstractQ) = rmul!(full!(A), Q) -rmul!(A::AbstractTriangular, Q::Adjoint{<:Any,<:AbstractQ}) = rmul!(full!(A), Q) - -function _qrmul(A, Q::AbstractQ) - TAQ = promote_type(eltype(A), eltype(Q)) - return rmul!(Matrix{TAQ}(A), convert(AbstractMatrix{TAQ}, Q)) -end -function _qrmul(A, adjQ::Adjoint{<:Any,<:AbstractQ}) - Q = adjQ.parent - TAQ = promote_type(eltype(A), eltype(Q)) - if size(A,2) == size(Q.factors, 1) - Anew = Matrix{TAQ}(A) - elseif size(A,2) == size(Q.factors,2) - Anew = [Matrix{TAQ}(A) zeros(TAQ, size(A, 1), size(Q.factors, 1) - size(Q.factors, 2))] - else - throw(DimensionMismatch("matrix A has dimensions $(size(A)) but matrix B has dimensions $(size(Q))")) - end - return rmul!(Anew, adjoint(convert(AbstractMatrix{TAQ}, Q))) -end +## Diagonal construction from UniformScaling +Diagonal{T}(s::UniformScaling, m::Integer) where {T} = Diagonal{T}(fill(T(s.λ), m)) +Diagonal(s::UniformScaling, m::Integer) = Diagonal{eltype(s)}(s, m) -*(A::AbstractTriangular, Q::AbstractQ) = _qrmul(A, Q) -*(A::AbstractTriangular, Q::Adjoint{<:Any,<:AbstractQ}) = _qrmul(A, Q) -*(A::BiTriSym, Q::AbstractQ) = _qrmul(A, Q) -*(A::BiTriSym, Q::Adjoint{<:Any,<:AbstractQ}) = _qrmul(A, Q) -*(A::Diagonal, Q::AbstractQ) = _qrmul(A, Q) -*(A::Diagonal, Q::Adjoint{<:Any,<:AbstractQ}) = _qrmul(A, Q) +Base.muladd(A::Union{Diagonal, UniformScaling}, B::Union{Diagonal, UniformScaling}, z::Union{Diagonal, UniformScaling}) = + Diagonal(_diag_or_value(A) .* _diag_or_value(B) .+ _diag_or_value(z)) -*(Q::AbstractQ, B::AbstractQ) = Q * (B * I) -*(Q::Adjoint{<:Any,<:AbstractQ}, B::AbstractQ) = Q * (B * I) -*(Q::AbstractQ, B::Adjoint{<:Any,<:AbstractQ}) = Q * (B * I) -*(Q::Adjoint{<:Any,<:AbstractQ}, B::Adjoint{<:Any,<:AbstractQ}) = Q * (B * I) +_diag_or_value(A::Diagonal) = A.diag +_diag_or_value(A::UniformScaling) = A.λ # fill[stored]! methods fillstored!(A::Diagonal, x) = (fill!(A.diag, x); A) diff --git a/stdlib/LinearAlgebra/src/uniformscaling.jl b/stdlib/LinearAlgebra/src/uniformscaling.jl index e691caf696e8c..f5c8bdd6f2e75 100644 --- a/stdlib/LinearAlgebra/src/uniformscaling.jl +++ b/stdlib/LinearAlgebra/src/uniformscaling.jl @@ -525,10 +525,6 @@ Array{T}(s::UniformScaling, m::Integer, n::Integer) where {T} = Matrix{T}(s, m, Array(s::UniformScaling, m::Integer, n::Integer) = Matrix(s, m, n) Array(s::UniformScaling, dims::Dims{2}) = Matrix(s, dims) -## Diagonal construction from UniformScaling -Diagonal{T}(s::UniformScaling, m::Integer) where {T} = Diagonal{T}(fill(T(s.λ), m)) -Diagonal(s::UniformScaling, m::Integer) = Diagonal{eltype(s)}(s, m) - dot(A::AbstractMatrix, J::UniformScaling) = dot(tr(A), J.λ) dot(J::UniformScaling, A::AbstractMatrix) = dot(J.λ, tr(A)) @@ -539,8 +535,3 @@ dot(x::AbstractVector, a::Union{Real,Complex}, y::AbstractVector) = a*dot(x, y) # muladd Base.muladd(A::UniformScaling, B::UniformScaling, z::UniformScaling) = UniformScaling(A.λ * B.λ + z.λ) -Base.muladd(A::Union{Diagonal, UniformScaling}, B::Union{Diagonal, UniformScaling}, z::Union{Diagonal, UniformScaling}) = - Diagonal(_diag_or_value(A) .* _diag_or_value(B) .+ _diag_or_value(z)) - -_diag_or_value(A::Diagonal) = A.diag -_diag_or_value(A::UniformScaling) = A.λ diff --git a/stdlib/LinearAlgebra/test/abstractq.jl b/stdlib/LinearAlgebra/test/abstractq.jl new file mode 100644 index 0000000000000..252a632fc97b9 --- /dev/null +++ b/stdlib/LinearAlgebra/test/abstractq.jl @@ -0,0 +1,84 @@ +# This file is a part of Julia. License is MIT: https://julialang.org/license + +module TestAbstractQ + +using Test +using LinearAlgebra +using LinearAlgebra: AbstractQ, AdjointQ +import LinearAlgebra: lmul!, rmul! +import Base: size, convert + +n = 5 + +@testset "custom AbstractQ type" begin + struct MyQ{T,S<:AbstractQ{T}} <: AbstractQ{T} + Q::S + end + MyQ{T}(Q::AbstractQ) where {T} = (P = convert(AbstractQ{T}, Q); MyQ{T,typeof(P)}(P)) + MyQ(Q::MyQ) = Q + + Base.size(Q::MyQ) = size(Q.Q) + LinearAlgebra.lmul!(Q::MyQ, B::AbstractVecOrMat) = lmul!(Q.Q, B) + LinearAlgebra.lmul!(adjQ::AdjointQ{<:Any,<:MyQ}, B::AbstractVecOrMat) = lmul!(parent(adjQ).Q', B) + LinearAlgebra.rmul!(A::AbstractMatrix, Q::MyQ) = rmul!(A, Q.Q) + LinearAlgebra.rmul!(A::AbstractMatrix, adjQ::AdjointQ{<:Any,<:MyQ}) = rmul!(A, parent(adjQ).Q') + Base.convert(::Type{AbstractQ{T}}, Q::MyQ) where {T} = MyQ{T}(Q.Q) + LinearAlgebra.det(Q::MyQ) = det(Q.Q) + + for T in (Float64, ComplexF64) + A = rand(T, n, n) + F = qr(A) + Q = MyQ(F.Q) + @test convert(AbstractQ{complex(T)}, Q) isa MyQ{complex(T)} + @test convert(AbstractQ{complex(T)}, Q') isa AdjointQ{<:complex(T),<:MyQ{complex(T)}} + @test Q*I ≈ Q.Q*I rtol=2eps(real(T)) + @test Q'*I ≈ Q.Q'*I rtol=2eps(real(T)) + @test I*Q ≈ Q.Q*I rtol=2eps(real(T)) + @test I*Q' ≈ I*Q.Q' rtol=2eps(real(T)) + @test abs(det(Q)) ≈ 1 + @test logabsdet(Q)[1] ≈ 0 atol=2eps(real(T)) + y = rand(T, n) + @test Q * y ≈ Q.Q * y ≈ Q' \ y ≈ ldiv!(Q', copy(y)) ≈ ldiv!(zero(y), Q', y) + @test Q'y ≈ Q.Q' * y ≈ Q \ y ≈ ldiv!(Q, copy(y)) ≈ ldiv!(zero(y), Q, y) + @test y'Q ≈ y'Q.Q ≈ y' / Q' + @test y'Q' ≈ y'Q.Q' ≈ y' / Q + y = Matrix(y') + @test y*Q ≈ y*Q.Q ≈ y / Q' ≈ rdiv!(copy(y), Q') + @test y*Q' ≈ y*Q.Q' ≈ y / Q ≈ rdiv!(copy(y), Q) + Y = rand(T, n, n); X = similar(Y) + for transQ in (identity, adjoint), transY in (identity, adjoint), Y in (Y, Y') + @test mul!(X, transQ(Q), transY(Y)) ≈ transQ(Q) * transY(Y) ≈ transQ(Q.Q) * transY(Y) + @test mul!(X, transY(Y), transQ(Q)) ≈ transY(Y) * transQ(Q) ≈ transY(Y) * transQ(Q.Q) + end + @test Matrix(Q) ≈ Q[:,:] ≈ copyto!(zeros(T, size(Q)), Q) ≈ Q.Q*I + @test Matrix(Q') ≈ (Q')[:,:] ≈ copyto!(zeros(T, size(Q)), Q') ≈ Q.Q'*I + @test Q[1,:] == Q.Q[1,:] + @test Q[:,1] == Q.Q[:,1] + @test Q[1,1] == Q.Q[1,1] + @test Q[:] == Q.Q[:] + @test Q[:,1:3] == Q.Q[:,1:3] == Matrix(Q)[:,1:3] + @test_throws BoundsError Q[0,1] + @test_throws BoundsError Q[n+1,1] + @test_throws BoundsError Q[1,0] + @test_throws BoundsError Q[1,n+1] + @test_throws BoundsError Q[:,1:n+1] + @test_throws BoundsError Q[:,0:n] + for perm in ((1, 2), (2, 1)) + P = PermutedDimsArray(zeros(T, size(Q)), perm) + @test copyto!(P, Q) ≈ Matrix(Q) + end + x = randn(T) + @test x * Q ≈ (x*I)*Q ≈ x * Q.Q + @test Q * x ≈ Q*(x*I) ≈ Q.Q * x + @test x * Q' ≈ (x*I)* Q' ≈ x * Q.Q' + @test Q' * x ≈ Q'*(x*I) ≈ Q.Q' * x + x = rand(T, 1) + Q = MyQ(qr(rand(T, 1, 1)).Q) + @test x * Q ≈ x * Q.Q + @test x * Q' ≈ x * Q.Q' + @test Q * x ≈ Q.Q * x + @test Q' * x ≈ Q.Q' * x + end +end + +end # module diff --git a/stdlib/LinearAlgebra/test/hessenberg.jl b/stdlib/LinearAlgebra/test/hessenberg.jl index fd1fefb97cab7..a9436b5ba8bdd 100644 --- a/stdlib/LinearAlgebra/test/hessenberg.jl +++ b/stdlib/LinearAlgebra/test/hessenberg.jl @@ -149,9 +149,11 @@ let n = 10 @test_throws ErrorException H.Z @test convert(Array, H) ≈ A @test (H.Q * H.H) * H.Q' ≈ A ≈ (Matrix(H.Q) * Matrix(H.H)) * Matrix(H.Q)' - @test (H.Q' *A) * H.Q ≈ H.H + @test (H.Q' * A) * H.Q ≈ H.H #getindex for HessenbergQ @test H.Q[1,1] ≈ Array(H.Q)[1,1] + @test det(H.Q) ≈ det(Matrix(H.Q)) + @test logabsdet(H.Q)[1] ≈ logabsdet(Matrix(H.Q))[1] atol=2n*eps(float(real(eltya))) # REPL show hessstring = sprint((t, s) -> show(t, "text/plain", s), H) diff --git a/stdlib/LinearAlgebra/test/lq.jl b/stdlib/LinearAlgebra/test/lq.jl index c340317a7cc23..8b4af6a0a5f8d 100644 --- a/stdlib/LinearAlgebra/test/lq.jl +++ b/stdlib/LinearAlgebra/test/lq.jl @@ -71,13 +71,11 @@ rectangularQ(Q::LinearAlgebra.LQPackedQ) = convert(Array, Q) @test lqa*x ≈ a*x rtol=3000ε @test (sq = size(q.factors, 2); *(Matrix{eltyb}(I, sq, sq), adjoint(q))*squareQ(q)) ≈ Matrix(I, n, n) rtol=5000ε if eltya != Int - @test Matrix{eltyb}(I, n, n)*q ≈ convert(AbstractMatrix{tab},q) + @test Matrix{eltyb}(I, n, n)*q ≈ Matrix(I, n, n) * convert(LinearAlgebra.AbstractQ{tab}, q) end @test q*x ≈ squareQ(q)*x rtol=100ε - @test transpose(q)*x ≈ transpose(squareQ(q))*x rtol=100ε @test q'*x ≈ squareQ(q)'*x rtol=100ε @test a*q ≈ a*squareQ(q) rtol=100ε - @test a*transpose(q) ≈ a*transpose(squareQ(q)) rtol=100ε @test a*q' ≈ a*squareQ(q)' rtol=100ε @test q*a'≈ squareQ(q)*a' rtol=100ε @test q'*a' ≈ squareQ(q)'*a' rtol=100ε @@ -89,7 +87,6 @@ rectangularQ(Q::LinearAlgebra.LQPackedQ) = convert(Array, Q) pad_a = vcat(I, a) pad_x = hcat(I, x) @test pad_a*q ≈ pad_a*squareQ(q) rtol=100ε - @test transpose(q)*pad_x ≈ transpose(squareQ(q))*pad_x rtol=100ε @test q'*pad_x ≈ squareQ(q)'*pad_x rtol=100ε end end @@ -193,12 +190,12 @@ end @testset for n in 1:3, m in 1:3 @testset "real" begin _, Q = lq(randn(n, m)) - @test det(Q) ≈ det(collect(Q)) + @test det(Q) ≈ det(Q*I) @test abs(det(Q)) ≈ 1 end @testset "complex" begin _, Q = lq(randn(ComplexF64, n, m)) - @test det(Q) ≈ det(collect(Q)) + @test det(Q) ≈ det(Q*I) @test abs(det(Q)) ≈ 1 end end @@ -217,11 +214,7 @@ L factor: 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 Q factor: -4×4 LinearAlgebra.LQPackedQ{Float64, Matrix{Float64}, Vector{Float64}}: - 1.0 0.0 0.0 0.0 - 0.0 1.0 0.0 0.0 - 0.0 0.0 1.0 0.0 - 0.0 0.0 0.0 1.0""" +4×4 LinearAlgebra.LQPackedQ{Float64, Matrix{Float64}, Vector{Float64}}""" end @testset "adjoint of LQ" begin diff --git a/stdlib/LinearAlgebra/test/qr.jl b/stdlib/LinearAlgebra/test/qr.jl index 9ccc42dfb3259..cb15d94d08dcb 100644 --- a/stdlib/LinearAlgebra/test/qr.jl +++ b/stdlib/LinearAlgebra/test/qr.jl @@ -21,8 +21,8 @@ breal = randn(n,2)/2 bimg = randn(n,2)/2 # helper functions to unambiguously recover explicit forms of an implicit QR Q -squareQ(Q::LinearAlgebra.AbstractQ) = (sq = size(Q.factors, 1); lmul!(Q, Matrix{eltype(Q)}(I, sq, sq))) -rectangularQ(Q::LinearAlgebra.AbstractQ) = convert(Array, Q) +squareQ(Q::LinearAlgebra.AbstractQ) = Q*I +rectangularQ(Q::LinearAlgebra.AbstractQ) = Matrix(Q) @testset for eltya in (Float32, Float64, ComplexF32, ComplexF64, BigFloat, Int) raw_a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal) @@ -62,7 +62,7 @@ rectangularQ(Q::LinearAlgebra.AbstractQ) = convert(Array, Q) sq = size(q.factors, 2) @test *(Matrix{eltyb}(I, sq, sq), adjoint(q)) * squareQ(q) ≈ Matrix(I, sq, sq) atol=5000ε if eltya != Int - @test Matrix{eltyb}(I, a_1, a_1)*q ≈ convert(AbstractMatrix{tab}, q) + @test Matrix{eltyb}(I, a_1, a_1)*q ≈ squareQ(convert(LinearAlgebra.AbstractQ{tab}, q)) ac = copy(a) @test qr!(a[:, 1:5])\b == qr!(view(ac, :, 1:5))\b end @@ -86,14 +86,14 @@ rectangularQ(Q::LinearAlgebra.AbstractQ) = convert(Array, Q) @test q*b[1:n1] ≈ rectangularQ(q)*b[1:n1] atol=100ε @test q*b ≈ squareQ(q)*b atol=100ε if eltya != Int - @test Array{eltya}(q) ≈ Matrix(q) + @test Array{eltya}(q) ≈ rectangularQ(q) end @test_throws DimensionMismatch q*b[1:n1 + 1] @test_throws DimensionMismatch b[1:n1 + 1]*q' sq = size(q.factors, 2) @test *(UpperTriangular(Matrix{eltyb}(I, sq, sq)), adjoint(q))*squareQ(q) ≈ Matrix(I, n1, a_1) atol=5000ε if eltya != Int - @test Matrix{eltyb}(I, a_1, a_1)*q ≈ convert(AbstractMatrix{tab},q) + @test Matrix{eltyb}(I, a_1, a_1)*q ≈ squareQ(convert(LinearAlgebra.AbstractQ{tab},q)) end # iterate q, r = qra @@ -123,7 +123,7 @@ rectangularQ(Q::LinearAlgebra.AbstractQ) = convert(Array, Q) @test_throws DimensionMismatch q*b[1:n1+1] @test_throws DimensionMismatch b[1:n1+1]*q' if eltya != Int - @test Matrix{eltyb}(I, n1, n1)*q ≈ convert(AbstractMatrix{tab},q) + @test Matrix{eltyb}(I, n1, n1)*q ≈ squareQ(convert(LinearAlgebra.AbstractQ{tab},q)) end # iterate q, r, p = qrpa @@ -149,7 +149,7 @@ rectangularQ(Q::LinearAlgebra.AbstractQ) = convert(Array, Q) sq = size(q.factors, 2) @test *(UpperTriangular(Matrix{eltyb}(I, sq, sq)), adjoint(q))*squareQ(q) ≈ Matrix(I, n1, a_1) atol=5000ε if eltya != Int - @test Matrix{eltyb}(I, a_1, a_1)*q ≈ convert(AbstractMatrix{tab},q) + @test Matrix{eltyb}(I, a_1, a_1)*q ≈ squareQ(convert(LinearAlgebra.AbstractQ{tab},q)) end qrstring = sprint((t, s) -> show(t, "text/plain", s), qrpa) rstring = sprint((t, s) -> show(t, "text/plain", s), r) @@ -205,6 +205,13 @@ rectangularQ(Q::LinearAlgebra.AbstractQ) = convert(Array, Q) @test mul!(c, b, q) ≈ b*q @test mul!(c, b, q') ≈ b*q' @test_throws DimensionMismatch mul!(Matrix{eltya}(I, n+1, n), q, b) + + b = similar(a[:,1]); rand!(b) + c = similar(a[:,1]) + d = similar(a[:,1]) + @test mul!(c, q, b) ≈ q*b + @test mul!(c, q', b) ≈ q'*b + @test_throws DimensionMismatch mul!(Vector{eltya}(undef, n+1), q, b) end end end @@ -228,7 +235,7 @@ end for T in (Tr, Complex{Tr}) v = convert(Vector{T}, vr) nv, nm = qr(v) - @test norm(nv - [-0.6 -0.8; -0.8 0.6], Inf) < eps(Tr) + @test norm(nv*Matrix(I, (2,2)) - [-0.6 -0.8; -0.8 0.6], Inf) < eps(Tr) @test nm == fill(-5.0, 1, 1) end end @@ -261,7 +268,7 @@ end @testset "Issue 24589. Promotion of rational matrices" begin A = rand(1//1:5//5, 4,3) - @test first(qr(A)) == first(qr(float(A))) + @test Matrix(first(qr(A))) == Matrix(first(qr(float(A)))) end @testset "Issue Test Factorization fallbacks for rectangular problems" begin @@ -303,7 +310,7 @@ end @testset for k in 0:min(n, m, 5) A = cat(Array(I(k)), randn(n - k, m - k); dims=(1, 2)) Q, = qr(A, pivot) - @test det(Q) ≈ det(collect(Q)) + @test det(Q) ≈ det(Q*Matrix(I, size(Q, 1), size(Q, 1))) @test abs(det(Q)) ≈ 1 end end @@ -311,7 +318,7 @@ end @testset for k in 0:min(n, m, 5) A = cat(Array(I(k)), randn(ComplexF64, n - k, m - k); dims=(1, 2)) Q, = qr(A, pivot) - @test det(Q) ≈ det(collect(Q)) + @test det(Q) ≈ det(Q*Matrix(I, size(Q, 1), size(Q, 1))) @test abs(det(Q)) ≈ 1 end end @@ -322,6 +329,7 @@ end for T in (Float64, ComplexF64) Q = qr(randn(T,5,5)).Q @test inv(Q) === Q' + @test inv(Q)' === inv(Q') === Q end end @@ -329,7 +337,7 @@ end for T in (Float32, Float64, ComplexF32, ComplexF64) Q1, R1 = qr(randn(T,5,5)) Q2, R2 = qr(Q1) - @test Q1 ≈ Q2 + @test Matrix(Q1) ≈ Matrix(Q2) @test R2 ≈ I end end @@ -362,13 +370,13 @@ end n = 5 Q, R = qr(randn(T,n,n)) Qmat = Matrix(Q) - dest1 = similar(Q) + dest1 = Matrix{T}(undef, size(Q)) copyto!(dest1, Q) @test dest1 ≈ Qmat - dest2 = PermutedDimsArray(similar(Q), (1, 2)) + dest2 = PermutedDimsArray(Matrix{T}(undef, size(Q)), (1, 2)) copyto!(dest2, Q) @test dest2 ≈ Qmat - dest3 = PermutedDimsArray(similar(Q), (2, 1)) + dest3 = PermutedDimsArray(Matrix{T}(undef, size(Q)), (2, 1)) copyto!(dest3, Q) @test dest3 ≈ Qmat end @@ -419,8 +427,8 @@ end A = qr(ones(3, 1)) B = I(3) C = B*A.Q' - @test C ≈ A.Q - @test A.Q' * B ≈ A.Q + @test C ≈ A.Q * Matrix(I, 3, 3) + @test A.Q' * B ≈ A.Q * Matrix(I, 3, 3) end @testset "convert between eltypes" begin diff --git a/stdlib/LinearAlgebra/test/testgroups b/stdlib/LinearAlgebra/test/testgroups index 2648016e453a8..8258573969cbb 100644 --- a/stdlib/LinearAlgebra/test/testgroups +++ b/stdlib/LinearAlgebra/test/testgroups @@ -25,4 +25,5 @@ bunchkaufman givens pinv factorization +abstractq ldlt diff --git a/test/bitarray.jl b/test/bitarray.jl index dd1d0d7d6c5a4..5d0bff62ab6e1 100644 --- a/test/bitarray.jl +++ b/test/bitarray.jl @@ -1652,7 +1652,7 @@ timesofar("cat") @test ((svdb1, svdb1A) = (svd(b1), svd(Array(b1))); svdb1.U == svdb1A.U && svdb1.S == svdb1A.S && svdb1.V == svdb1A.V) @test ((qrb1, qrb1A) = (qr(b1), qr(Array(b1))); - qrb1.Q == qrb1A.Q && qrb1.R == qrb1A.R) + Matrix(qrb1.Q) == Matrix(qrb1A.Q) && qrb1.R == qrb1A.R) b1 = bitrand(v1) @check_bit_operation diagm(0 => b1) BitMatrix