diff --git a/base/linalg/qr.jl b/base/linalg/qr.jl index 6cdd2d0dd3ede..36c9e14d13e81 100644 --- a/base/linalg/qr.jl +++ b/base/linalg/qr.jl @@ -422,6 +422,8 @@ function getindex(A::QRPivoted{T}, d::Symbol) where T end end +abstract type AbstractQ{T} <: AbstractMatrix{T} end + # Type-stable interface to get Q getq(A::QRCompactWY) = QRCompactWYQ(A.factors,A.T) getq(A::Union{QR, QRPivoted}) = QRPackedQ(A.factors,A.τ) @@ -432,7 +434,7 @@ getq(A::Union{QR, QRPivoted}) = QRPackedQ(A.factors,A.τ) The orthogonal/unitary ``Q`` matrix of a QR factorization stored in [`QR`](@ref) or [`QRPivoted`](@ref) format. """ -struct QRPackedQ{T,S<:AbstractMatrix} <: AbstractMatrix{T} +struct QRPackedQ{T,S<:AbstractMatrix} <: AbstractQ{T} factors::S τ::Vector{T} QRPackedQ{T,S}(factors::AbstractMatrix{T}, τ::Vector{T}) where {T,S<:AbstractMatrix} = new(factors, τ) @@ -445,7 +447,7 @@ QRPackedQ(factors::AbstractMatrix{T}, τ::Vector{T}) where {T} = QRPackedQ{T,typ The orthogonal/unitary ``Q`` matrix of a QR factorization stored in [`QRCompactWY`](@ref) format. """ -struct QRCompactWYQ{S, M<:AbstractMatrix} <: AbstractMatrix{S} +struct QRCompactWYQ{S, M<:AbstractMatrix} <: AbstractQ{S} factors::M T::Matrix{S} QRCompactWYQ{S,M}(factors::AbstractMatrix{S}, T::Matrix{S}) where {S,M<:AbstractMatrix} = new(factors, T) @@ -458,11 +460,11 @@ convert(::Type{AbstractMatrix{T}}, Q::QRPackedQ) where {T} = convert(QRPackedQ{T convert(::Type{QRCompactWYQ{S}}, Q::QRCompactWYQ) where {S} = QRCompactWYQ(convert(AbstractMatrix{S}, Q.factors), convert(AbstractMatrix{S}, Q.T)) convert(::Type{AbstractMatrix{S}}, Q::QRCompactWYQ{S}) where {S} = Q convert(::Type{AbstractMatrix{S}}, Q::QRCompactWYQ) where {S} = convert(QRCompactWYQ{S}, Q) -convert(::Type{Matrix}, A::Union{QRPackedQ{T},QRCompactWYQ{T}}) where {T} = A_mul_B!(A, eye(T, size(A.factors, 1), minimum(size(A.factors)))) -convert(::Type{Array}, A::Union{QRPackedQ,QRCompactWYQ}) = convert(Matrix, A) +convert(::Type{Matrix}, A::AbstractQ{T}) where {T} = A_mul_B!(A, eye(T, size(A.factors, 1), minimum(size(A.factors)))) +convert(::Type{Array}, A::AbstractQ) = convert(Matrix, A) """ - full(A::Union{QRPackedQ,QRCompactWYQ}; thin::Bool=true) -> Matrix + full(A::AbstractQ; thin::Bool=true) -> Matrix Converts an orthogonal or unitary matrix stored as a `QRCompactWYQ` object, i.e. in the compact WY format [^Bischof1987], or in the `QRPackedQ` format, to a dense matrix. @@ -472,7 +474,7 @@ rows of `R` in the QR factorization that are zero. The resulting matrix is the ` QR factorization (sometimes called the reduced QR factorization). If `false`, returns a `Q` that spans all rows of `R` in its corresponding QR factorization. """ -function full(A::Union{QRPackedQ{T},QRCompactWYQ{T}}; thin::Bool = true) where T +function full(A::AbstractQ{T}; thin::Bool = true) where T if thin convert(Array, A) else @@ -482,11 +484,11 @@ end size(A::Union{QR,QRCompactWY,QRPivoted}, dim::Integer) = size(A.factors, dim) size(A::Union{QR,QRCompactWY,QRPivoted}) = size(A.factors) -size(A::Union{QRPackedQ,QRCompactWYQ}, dim::Integer) = 0 < dim ? (dim <= 2 ? size(A.factors, 1) : 1) : throw(BoundsError()) -size(A::Union{QRPackedQ,QRCompactWYQ}) = size(A, 1), size(A, 2) +size(A::AbstractQ, dim::Integer) = 0 < dim ? (dim <= 2 ? size(A.factors, 1) : 1) : throw(BoundsError()) +size(A::AbstractQ) = size(A, 1), size(A, 2) -function getindex(A::Union{QRPackedQ,QRCompactWYQ}, i::Integer, j::Integer) +function getindex(A::AbstractQ, i::Integer, j::Integer) x = zeros(eltype(A), size(A, 1)) x[i] = 1 y = zeros(eltype(A), size(A, 2)) @@ -496,8 +498,10 @@ end ## Multiplication by Q ### QB -A_mul_B!(A::QRCompactWYQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = LAPACK.gemqrt!('L','N',A.factors,A.T,B) -A_mul_B!(A::QRPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = LAPACK.ormqr!('L','N',A.factors,A.τ,B) +A_mul_B!(A::QRCompactWYQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasFloat, S<:StridedMatrix} = + LAPACK.gemqrt!('L','N',A.factors,A.T,B) +A_mul_B!(A::QRPackedQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasFloat, S<:StridedMatrix} = + LAPACK.ormqr!('L','N',A.factors,A.τ,B) function A_mul_B!(A::QRPackedQ, B::AbstractVecOrMat) mA, nA = size(A.factors) mB, nB = size(B,1), size(B,2) @@ -523,7 +527,7 @@ function A_mul_B!(A::QRPackedQ, B::AbstractVecOrMat) B end -function (*)(A::Union{QRPackedQ,QRCompactWYQ}, b::StridedVector) +function (*)(A::AbstractQ, b::StridedVector) TAb = promote_type(eltype(A), eltype(b)) Anew = convert(AbstractMatrix{TAb}, A) if size(A.factors, 1) == length(b) @@ -535,7 +539,7 @@ function (*)(A::Union{QRPackedQ,QRCompactWYQ}, b::StridedVector) end A_mul_B!(Anew, bnew) end -function (*)(A::Union{QRPackedQ,QRCompactWYQ}, B::StridedMatrix) +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) @@ -549,10 +553,14 @@ function (*)(A::Union{QRPackedQ,QRCompactWYQ}, B::StridedMatrix) end ### QcB -Ac_mul_B!(A::QRCompactWYQ{T}, B::StridedVecOrMat{T}) where {T<:BlasReal} = LAPACK.gemqrt!('L','T',A.factors,A.T,B) -Ac_mul_B!(A::QRCompactWYQ{T}, B::StridedVecOrMat{T}) where {T<:BlasComplex} = LAPACK.gemqrt!('L','C',A.factors,A.T,B) -Ac_mul_B!(A::QRPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasReal} = LAPACK.ormqr!('L','T',A.factors,A.τ,B) -Ac_mul_B!(A::QRPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasComplex} = LAPACK.ormqr!('L','C',A.factors,A.τ,B) +Ac_mul_B!(A::QRCompactWYQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasReal,S<:StridedMatrix} = + LAPACK.gemqrt!('L','T',A.factors,A.T,B) +Ac_mul_B!(A::QRCompactWYQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasComplex,S<:StridedMatrix} = + LAPACK.gemqrt!('L','C',A.factors,A.T,B) +Ac_mul_B!(A::QRPackedQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasReal,S<:StridedMatrix} = + LAPACK.ormqr!('L','T',A.factors,A.τ,B) +Ac_mul_B!(A::QRPackedQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasComplex,S<:StridedMatrix} = + LAPACK.ormqr!('L','C',A.factors,A.τ,B) function Ac_mul_B!(A::QRPackedQ, B::AbstractVecOrMat) mA, nA = size(A.factors) mB, nB = size(B,1), size(B,2) @@ -577,7 +585,7 @@ function Ac_mul_B!(A::QRPackedQ, B::AbstractVecOrMat) end B end -function Ac_mul_B(Q::Union{QRPackedQ,QRCompactWYQ}, B::StridedVecOrMat) +function Ac_mul_B(Q::AbstractQ, B::StridedVecOrMat) TQB = promote_type(eltype(Q), eltype(B)) return Ac_mul_B!(convert(AbstractMatrix{TQB}, Q), copy_oftype(B, TQB)) end @@ -586,7 +594,7 @@ end for (f1, f2) in ((:A_mul_Bc, :A_mul_B!), (:Ac_mul_Bc, :Ac_mul_B!)) @eval begin - function ($f1)(Q::Union{QRPackedQ,QRCompactWYQ}, B::StridedVecOrMat) + function ($f1)(Q::AbstractQ, B::StridedVecOrMat) TQB = promote_type(eltype(Q), eltype(B)) Bc = similar(B, TQB, (size(B, 2), size(B, 1))) ctranspose!(Bc, B) @@ -596,8 +604,10 @@ for (f1, f2) in ((:A_mul_Bc, :A_mul_B!), end ### AQ -A_mul_B!(A::StridedVecOrMat{T}, B::QRCompactWYQ{T}) where {T<:BlasFloat} = LAPACK.gemqrt!('R','N', B.factors, B.T, A) -A_mul_B!(A::StridedVecOrMat{T}, B::QRPackedQ{T}) where {T<:BlasFloat} = LAPACK.ormqr!('R', 'N', B.factors, B.τ, A) +A_mul_B!(A::StridedVecOrMat{T}, B::QRCompactWYQ{T,S}) where {T<:BlasFloat,S<:StridedMatrix} = + LAPACK.gemqrt!('R','N', B.factors, B.T, A) +A_mul_B!(A::StridedVecOrMat{T}, B::QRPackedQ{T,S}) where {T<:BlasFloat,S<:StridedMatrix} = + LAPACK.ormqr!('R', 'N', B.factors, B.τ, A) function A_mul_B!(A::StridedMatrix,Q::QRPackedQ) mQ, nQ = size(Q.factors) mA, nA = size(A,1), size(A,2) @@ -623,7 +633,7 @@ function A_mul_B!(A::StridedMatrix,Q::QRPackedQ) A end -function (*)(A::StridedMatrix, Q::Union{QRPackedQ,QRCompactWYQ}) +function (*)(A::StridedMatrix, Q::AbstractQ) TAQ = promote_type(eltype(A), eltype(Q)) return A_mul_B!(copy_oftype(A, TAQ), convert(AbstractMatrix{TAQ}, Q)) end @@ -633,7 +643,7 @@ A_mul_Bc!(A::StridedVecOrMat{T}, B::QRCompactWYQ{T}) where {T<:BlasReal} = LAPAC A_mul_Bc!(A::StridedVecOrMat{T}, B::QRCompactWYQ{T}) where {T<:BlasComplex} = LAPACK.gemqrt!('R','C',B.factors,B.T,A) A_mul_Bc!(A::StridedVecOrMat{T}, B::QRPackedQ{T}) where {T<:BlasReal} = LAPACK.ormqr!('R','T',B.factors,B.τ,A) A_mul_Bc!(A::StridedVecOrMat{T}, B::QRPackedQ{T}) where {T<:BlasComplex} = LAPACK.ormqr!('R','C',B.factors,B.τ,A) -function A_mul_Bc!(A::AbstractMatrix,Q::QRPackedQ) +function A_mul_Bc!(A::StridedMatrix,Q::QRPackedQ) mQ, nQ = size(Q.factors) mA, nA = size(A,1), size(A,2) if nA != mQ @@ -657,7 +667,7 @@ function A_mul_Bc!(A::AbstractMatrix,Q::QRPackedQ) end A end -function A_mul_Bc(A::AbstractMatrix, B::Union{QRCompactWYQ,QRPackedQ}) +function A_mul_Bc(A::StridedMatrix, B::AbstractQ) TAB = promote_type(eltype(A),eltype(B)) BB = convert(AbstractMatrix{TAB}, B) if size(A,2) == size(B.factors, 1) @@ -670,14 +680,14 @@ function A_mul_Bc(A::AbstractMatrix, B::Union{QRCompactWYQ,QRPackedQ}) throw(DimensionMismatch("matrix A has dimensions $(size(A)) but matrix B has dimensions $(size(B))")) end end -@inline A_mul_Bc(rowvec::RowVector, B::Union{LinAlg.QRCompactWYQ,LinAlg.QRPackedQ}) = ctranspose(B*ctranspose(rowvec)) +@inline A_mul_Bc(rowvec::RowVector, B::AbstractQ) = ctranspose(B*ctranspose(rowvec)) ### AcQ/AcQc for (f1, f2) in ((:Ac_mul_B, :A_mul_B!), (:Ac_mul_Bc, :A_mul_Bc!)) @eval begin - function ($f1)(A::StridedVecOrMat, Q::Union{QRPackedQ,QRCompactWYQ}) + function ($f1)(A::StridedVecOrMat, Q::AbstractQ) TAQ = promote_type(eltype(A), eltype(Q)) Ac = similar(A, TAQ, (size(A, 2), size(A, 1))) ctranspose!(Ac, A) diff --git a/base/sparse/cholmod.jl b/base/sparse/cholmod.jl index f049dea87c553..5954638bbf22e 100644 --- a/base/sparse/cholmod.jl +++ b/base/sparse/cholmod.jl @@ -248,6 +248,8 @@ mutable struct Sparse{Tv<:VTypes} <: AbstractSparseMatrix{Tv,SuiteSparse_long} end Sparse(p::Ptr{C_Sparse{Tv}}) where {Tv<:VTypes} = Sparse{Tv}(p) +Base.unsafe_convert(::Type{Ptr{Tv}}, A::Sparse{Tv}) where {Tv} = A.p + # Factor if build_version >= v"2.1.0" # CHOLMOD version 2.1.0 or later diff --git a/base/sparse/linalg.jl b/base/sparse/linalg.jl index d0116fb5bcdb8..eefd80f9d09b6 100644 --- a/base/sparse/linalg.jl +++ b/base/sparse/linalg.jl @@ -194,7 +194,7 @@ function spmatmul(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti}; end ## solvers -function fwdTriSolve!(A::SparseMatrixCSC, B::AbstractVecOrMat) +function fwdTriSolve!(A::SparseMatrixCSCUnion, B::AbstractVecOrMat) # forward substitution for CSC matrices nrowB, ncolB = size(B, 1), size(B, 2) ncol = LinAlg.checksquare(A) @@ -202,9 +202,9 @@ function fwdTriSolve!(A::SparseMatrixCSC, B::AbstractVecOrMat) throw(DimensionMismatch("A is $(ncol) columns and B has $(nrowB) rows")) end - aa = A.nzval - ja = A.rowval - ia = A.colptr + aa = getnzval(A) + ja = getrowval(A) + ia = getcolptr(A) joff = 0 for k = 1:ncolB @@ -239,7 +239,7 @@ function fwdTriSolve!(A::SparseMatrixCSC, B::AbstractVecOrMat) B end -function bwdTriSolve!(A::SparseMatrixCSC, B::AbstractVecOrMat) +function bwdTriSolve!(A::SparseMatrixCSCUnion, B::AbstractVecOrMat) # backward substitution for CSC matrices nrowB, ncolB = size(B, 1), size(B, 2) ncol = LinAlg.checksquare(A) @@ -247,9 +247,9 @@ function bwdTriSolve!(A::SparseMatrixCSC, B::AbstractVecOrMat) throw(DimensionMismatch("A is $(ncol) columns and B has $(nrowB) rows")) end - aa = A.nzval - ja = A.rowval - ia = A.colptr + aa = getnzval(A) + ja = getrowval(A) + ia = getcolptr(A) joff = 0 for k = 1:ncolB @@ -284,11 +284,11 @@ function bwdTriSolve!(A::SparseMatrixCSC, B::AbstractVecOrMat) B end -A_ldiv_B!(L::LowerTriangular{T,<:SparseMatrixCSC{T}}, B::StridedVecOrMat) where {T} = fwdTriSolve!(L.data, B) -A_ldiv_B!(U::UpperTriangular{T,<:SparseMatrixCSC{T}}, B::StridedVecOrMat) where {T} = bwdTriSolve!(U.data, B) +A_ldiv_B!(L::LowerTriangular{T,<:SparseMatrixCSCUnion{T}}, B::StridedVecOrMat) where {T} = fwdTriSolve!(L.data, B) +A_ldiv_B!(U::UpperTriangular{T,<:SparseMatrixCSCUnion{T}}, B::StridedVecOrMat) where {T} = bwdTriSolve!(U.data, B) -(\)(L::LowerTriangular{T,<:SparseMatrixCSC{T}}, B::SparseMatrixCSC) where {T} = A_ldiv_B!(L, Array(B)) -(\)(U::UpperTriangular{T,<:SparseMatrixCSC{T}}, B::SparseMatrixCSC) where {T} = A_ldiv_B!(U, Array(B)) +(\)(L::LowerTriangular{T,<:SparseMatrixCSCUnion{T}}, B::SparseMatrixCSC) where {T} = A_ldiv_B!(L, Array(B)) +(\)(U::UpperTriangular{T,<:SparseMatrixCSCUnion{T}}, B::SparseMatrixCSC) where {T} = A_ldiv_B!(U, Array(B)) function A_rdiv_B!(A::SparseMatrixCSC{T}, D::Diagonal{T}) where T dd = D.diag diff --git a/base/sparse/sparsematrix.jl b/base/sparse/sparsematrix.jl index 709016fa9678c..c474c7006f1f2 100644 --- a/base/sparse/sparsematrix.jl +++ b/base/sparse/sparsematrix.jl @@ -35,6 +35,24 @@ end size(S::SparseMatrixCSC) = (S.m, S.n) +# Define an alias for views of a SparseMatrixCSC which include all rows and a unit range of the columns. +# Also define a union of SparseMatrixCSC and this view since many methods can be defined efficiently for +# this union by extracting the fields via the get function: getcolptr, getrowval, and getnzval. The key +# insight is that getcolptr on a SparseMatrixCSCView returns an offset view of the colptr of the +# underlying SparseMatrixCSC +const SparseMatrixCSCView{Tv,Ti} = + SubArray{Tv,2,SparseMatrixCSC{Tv,Ti}, + Tuple{Base.Slice{Base.OneTo{Int}},I}} where {I<:AbstractUnitRange} +const SparseMatrixCSCUnion{Tv,Ti} = Union{SparseMatrixCSC{Tv,Ti}, SparseMatrixCSCView{Tv,Ti}} + +getcolptr(S::SparseMatrixCSC) = S.colptr +getcolptr(S::SparseMatrixCSCView) = view(S.parent.colptr, first(indices(S, 2)):(last(indices(S, 2)) + 1)) +getrowval(S::SparseMatrixCSC) = S.rowval +getrowval(S::SparseMatrixCSCView) = S.parent.rowval +getnzval( S::SparseMatrixCSC) = S.nzval +getnzval( S::SparseMatrixCSCView) = S.parent.nzval + + """ nnz(A) diff --git a/base/sparse/sparsevector.jl b/base/sparse/sparsevector.jl index 73765f5cfb9d4..784f462cc8dd1 100644 --- a/base/sparse/sparsevector.jl +++ b/base/sparse/sparsevector.jl @@ -30,6 +30,11 @@ end SparseVector(n::Integer, nzind::Vector{Ti}, nzval::Vector{Tv}) where {Tv,Ti} = SparseVector{Tv,Ti}(n, nzind, nzval) +# Define an alias for a view of a whole column of a SparseMatrixCSC. Many methods can be written for the +# union of such a view and a SparseVector so we define an alias for such a union as well +const SparseColumnView{T} = SubArray{T,1,<:SparseMatrixCSC,Tuple{Base.Slice{Base.OneTo{Int}},Int},false} +const SparseVectorUnion{T} = Union{SparseVector{T}, SparseColumnView{T}} + ### Basic properties length(x::SparseVector) = x.n @@ -37,8 +42,22 @@ size(x::SparseVector) = (x.n,) nnz(x::SparseVector) = length(x.nzval) countnz(x::SparseVector) = countnz(x.nzval) count(x::SparseVector) = count(x.nzval) + nonzeros(x::SparseVector) = x.nzval +function nonzeros(x::SparseColumnView) + rowidx, colidx = parentindexes(x) + A = parent(x) + @inbounds y = view(A.nzval, nzrange(A, colidx)) + return y +end + nonzeroinds(x::SparseVector) = x.nzind +function nonzeroinds(x::SparseColumnView) + rowidx, colidx = parentindexes(x) + A = parent(x) + @inbounds y = view(A.rowval, nzrange(A, colidx)) + return y +end similar(x::SparseVector, Tv::Type=eltype(x)) = SparseVector(x.n, copy(x.nzind), Vector{Tv}(length(x.nzval))) function similar(x::SparseVector, ::Type{Tv}, ::Type{Ti}) where {Tv,Ti} @@ -1396,7 +1415,7 @@ for f in [:sum, :maximum, :minimum], op in [:abs, :abs2] end end -vecnorm(x::AbstractSparseVector, p::Real=2) = vecnorm(nonzeros(x), p) +vecnorm(x::SparseVectorUnion, p::Real=2) = vecnorm(nonzeros(x), p) ### linalg.jl @@ -1409,7 +1428,7 @@ vecnorm(x::AbstractSparseVector, p::Real=2) = vecnorm(nonzeros(x), p) # axpy -function LinAlg.axpy!(a::Number, x::AbstractSparseVector, y::StridedVector) +function LinAlg.axpy!(a::Number, x::SparseVectorUnion, y::StridedVector) length(x) == length(y) || throw(DimensionMismatch()) nzind = nonzeroinds(x) nzval = nonzeros(x) @@ -1454,8 +1473,7 @@ broadcast(::typeof(*), a::Number, x::AbstractSparseVector) = a * x broadcast(::typeof(/), x::AbstractSparseVector, a::Number) = x / a # dot - -function dot(x::StridedVector{Tx}, y::AbstractSparseVector{Ty}) where {Tx<:Number,Ty<:Number} +function dot(x::StridedVector{Tx}, y::SparseVectorUnion{Ty}) where {Tx<:Number,Ty<:Number} n = length(x) length(y) == n || throw(DimensionMismatch()) nzind = nonzeroinds(y) @@ -1467,13 +1485,13 @@ function dot(x::StridedVector{Tx}, y::AbstractSparseVector{Ty}) where {Tx<:Numbe return s end -function dot(x::AbstractSparseVector{Tx}, y::AbstractVector{Ty}) where {Tx<:Number,Ty<:Number} +function dot(x::SparseVectorUnion{Tx}, y::AbstractVector{Ty}) where {Tx<:Number,Ty<:Number} n = length(y) length(x) == n || throw(DimensionMismatch()) nzind = nonzeroinds(x) nzval = nonzeros(x) s = zero(Tx) * zero(Ty) - for i = 1:length(nzind) + @inbounds for i = 1:length(nzind) s += conj(nzval[i]) * y[nzind[i]] end return s @@ -1500,7 +1518,7 @@ function _spdot(f::Function, s end -function dot(x::AbstractSparseVector{<:Number}, y::AbstractSparseVector{<:Number}) +function dot(x::SparseVectorUnion{<:Number}, y::SparseVectorUnion{<:Number}) x === y && return sum(abs2, x) n = length(x) length(y) == n || throw(DimensionMismatch()) @@ -1518,6 +1536,19 @@ end ### BLAS-2 / dense A * sparse x -> dense y +# lowrankupdate (BLAS.ger! like) +function LinAlg.lowrankupdate!(A::StridedMatrix, x::StridedVector, y::SparseVectorUnion, α::Number = 1) + nzi = nonzeroinds(y) + nzv = nonzeros(y) + @inbounds for (j,v) in zip(nzi,nzv) + αv = α*v' + for i in indices(x, 1) + A[i,j] += x[i]*αv + end + end + return A +end + # A_mul_B function (*)(A::StridedMatrix{Ta}, x::AbstractSparseVector{Tx}) where {Ta,Tx} diff --git a/base/sparse/spqr.jl b/base/sparse/spqr.jl index a507c8cdcbf22..b763a6d10cb87 100644 --- a/base/sparse/spqr.jl +++ b/base/sparse/spqr.jl @@ -19,135 +19,322 @@ const ORDERING_BESTAMD = Int32(9) # try COLAMD and AMD; pick best# # Let [m n] = size of the matrix after pruning singletons. The default # ordering strategy is to use COLAMD if m <= 2*n. Otherwise, AMD(A'A) is # tried. If there is a high fill-in with AMD then try METIS(A'A) and take -# the best of AMD and METIS. METIS is not tried if it isn't installed. +# the best of AMD and METIS. METIS is not tried if it isn't installed. -# tol options -const DEFAULT_TOL = Int32(-2) # if tol <= -2, the default tol is used -const NO_TOL = Int32(-1) # if -2 < tol < 0, then no tol is used - -# for qmult, method can be 0,1,2,3: -const QTX = Int32(0) -const QX = Int32(1) -const XQT = Int32(2) -const XQ = Int32(3) +using ..SparseArrays: SparseMatrixCSC +using ..SparseArrays.CHOLMOD +using ..SparseArrays.CHOLMOD: change_stype!, free! -# system can be 0,1,2,3: Given Q*R=A*E from SuiteSparseQR_factorize: -const RX_EQUALS_B = Int32(0) # solve R*X=B or X = R\B -const RETX_EQUALS_B = Int32(1) # solve R*E'*X=B or X = E*(R\B) -const RTX_EQUALS_B = Int32(2) # solve R'*X=B or X = R'\B -const RTX_EQUALS_ETB = Int32(3) # solve R'*X=E'*B or X = R'\(E'*B) +function _qr!(ordering::Integer, tol::Real, econ::Integer, getCTX::Integer, + A::Sparse{Tv}, + Bsparse::Union{Sparse{Tv} , Ptr{Void}} = C_NULL, + Bdense::Union{Dense{Tv} , Ptr{Void}} = C_NULL, + Zsparse::Union{Ref{Ptr{CHOLMOD.C_Sparse{Tv}}} , Ptr{Void}} = C_NULL, + Zdense::Union{Ref{Ptr{CHOLMOD.C_Dense{Tv}}} , Ptr{Void}} = C_NULL, + R::Union{Ref{Ptr{CHOLMOD.C_Sparse{Tv}}} , Ptr{Void}} = C_NULL, + E::Union{Ref{Ptr{CHOLMOD.SuiteSparse_long}} , Ptr{Void}} = C_NULL, + H::Union{Ref{Ptr{CHOLMOD.C_Sparse{Tv}}} , Ptr{Void}} = C_NULL, + HPinv::Union{Ref{Ptr{CHOLMOD.SuiteSparse_long}}, Ptr{Void}} = C_NULL, + HTau::Union{Ref{Ptr{CHOLMOD.C_Dense{Tv}}} , Ptr{Void}} = C_NULL) where {Tv<:CHOLMOD.VTypes} + AA = unsafe_load(pointer(A)) + m, n = AA.nrow, AA.ncol + rnk = ccall((:SuiteSparseQR_C, :libspqr), CHOLMOD.SuiteSparse_long, + (Cint, Cdouble, CHOLMOD.SuiteSparse_long, Cint, + Ptr{CHOLMOD.C_Sparse{Tv}}, Ptr{CHOLMOD.C_Sparse{Tv}}, Ptr{CHOLMOD.C_Dense{Tv}}, + Ptr{Ptr{CHOLMOD.C_Sparse{Tv}}}, Ptr{Ptr{CHOLMOD.C_Dense{Tv}}}, Ptr{Ptr{CHOLMOD.C_Sparse{Tv}}}, + Ptr{Ptr{CHOLMOD.SuiteSparse_long}}, Ptr{Ptr{CHOLMOD.C_Sparse{Tv}}}, Ptr{Ptr{CHOLMOD.SuiteSparse_long}}, + Ptr{Ptr{CHOLMOD.C_Dense{Tv}}}, Ptr{Void}), + ordering, # all, except 3:given treated as 0:fixed + tol, # columns with 2-norm <= tol treated as 0 + econ, # e = max(min(m,econ),rank(A)) + getCTX, # 0: Z=C (e-by-k), 1: Z=C', 2: Z=X (e-by-k) + A, # m-by-n sparse matrix to factorize + Bsparse, # sparse m-by-k B + Bdense, # dense m-by-k B + # /* outputs: */ + Zsparse, # sparse Z + Zdense, # dense Z + R, # e-by-n sparse matrix */ + E, # size n column perm, NULL if identity */ + H, # m-by-nh Householder vectors + HPinv, # size m row permutation + HTau, # 1-by-nh Householder coefficients + CHOLMOD.common()) # /* workspace and parameters */ -using ..SparseArrays: SparseMatrixCSC -using ..SparseArrays.CHOLMOD: C_Dense, C_Sparse, Dense, ITypes, Sparse, SuiteSparseStruct, VTypes, common + if rnk < 0 + error("Sparse QR factorization failed") + end -import Base: size -import Base.LinAlg: qrfact -import ..SparseArrays.CHOLMOD: convert, free! + e = E[] + if e == C_NULL + _E = Vector{CHOLMOD.SuiteSparse_long}() + else + _E = Vector{CHOLMOD.SuiteSparse_long}(n) + for i in 1:n + @inbounds _E[i] = unsafe_load(e, i) + 1 + end + # Free memory allocated by SPQR. This call will make sure that the + # correct deallocator function is called and that the memory count in + # the common struct is updated + ccall((:cholmod_l_free, :libcholmod), Void, + (Csize_t, Cint, Ptr{CHOLMOD.SuiteSparse_long}, Ptr{Void}), + n, sizeof(CHOLMOD.SuiteSparse_long), e, CHOLMOD.common()) + end + hpinv = HPinv[] + if hpinv == C_NULL + _HPinv = Vector{CHOLMOD.SuiteSparse_long}() + else + _HPinv = Vector{CHOLMOD.SuiteSparse_long}(m) + for i in 1:m + @inbounds _HPinv[i] = unsafe_load(hpinv, i) + 1 + end + # Free memory allocated by SPQR. This call will make sure that the + # correct deallocator function is called and that the memory count in + # the common struct is updated + ccall((:cholmod_l_free, :libcholmod), Void, + (Csize_t, Cint, Ptr{CHOLMOD.SuiteSparse_long}, Ptr{Void}), + m, sizeof(CHOLMOD.SuiteSparse_long), hpinv, CHOLMOD.common()) + end -struct C_Factorization{Tv<:VTypes} <: SuiteSparseStruct - xtype::Cint - factors::Ptr{Tv} + return rnk, _E, _HPinv end -mutable struct Factorization{Tv<:VTypes} <: Base.LinAlg.Factorization{Tv} - m::Int - n::Int - p::Ptr{C_Factorization{Tv}} - function Factorization{Tv}(m::Integer, n::Integer, p::Ptr{C_Factorization{Tv}}) where Tv<:VTypes - if p == C_NULL - throw(ArgumentError("factorization failed for unknown reasons. Please submit a bug report.")) - end - new(m, n, p) - end +# Struct for storing sparse QR from SPQR such that +# A[invperm(rpivinv), cpiv] = (I - factors[:,1]*τ[1]*factors[:,1]')*...*(I - factors[:,k]*τ[k]*factors[:,k]')*R +# with k = size(factors, 2). +struct QRSparse{Tv,Ti} <: LinAlg.Factorization{Tv} + factors::SparseMatrixCSC{Tv,Ti} + τ::Vector{Tv} + R::SparseMatrixCSC{Tv,Ti} + cpiv::Vector{Ti} + rpivinv::Vector{Ti} end -Factorization(m::Integer, n::Integer, p::Ptr{C_Factorization{Tv}}) where Tv<:VTypes = Factorization{Tv}(m, n, p) -size(F::Factorization) = (F.m, F.n) -function size(F::Factorization, i::Integer) - if i < 1 - throw(ArgumentError("dimension must be positive")) - end +Base.size(F::QRSparse) = (size(F.factors, 1), size(F.R, 2)) +function Base.size(F::QRSparse, i::Integer) if i == 1 - return F.m + return size(F.factors, 1) elseif i == 2 - return F.n + return size(F.R, 2) + elseif i > 2 + return 1 + else + throw(ArgumentError("second argument must be positive")) end - return 1 end -function free!(F::Factorization{Tv}) where Tv<:VTypes - ccall((:SuiteSparseQR_C_free, :libspqr), Cint, - (Ptr{Ptr{C_Factorization{Tv}}}, Ptr{Void}), - &F.p, common()) == 1 +struct QRSparseQ{Tv<:CHOLMOD.VTypes,Ti<:Integer} <: LinAlg.AbstractQ{Tv} + factors::SparseMatrixCSC{Tv,Ti} + τ::Vector{Tv} end -function backslash(ordering::Integer, tol::Real, A::Sparse{Tv}, B::Dense{Tv}) where Tv<:VTypes - m, n = size(A) - if m != size(B, 1) - throw(DimensionMismatch("left hand side and right hand side must have same number of rows")) - end - d = Dense(ccall((:SuiteSparseQR_C_backslash, :libspqr), Ptr{C_Dense{Tv}}, - (Cint, Cdouble, Ptr{C_Sparse{Tv}}, Ptr{C_Dense{Tv}}, Ptr{Void}), - ordering, tol, get(A.p), get(B.p), common())) - finalizer(d, free!) - d +Base.size(Q::QRSparseQ) = (size(Q.factors, 1), size(Q.factors, 1)) + +# From SPQR manual p. 6 +_default_tol(A::SparseMatrixCSC) = + 20*sum(size(A))*eps(real(eltype(A)))*maximum(norm(view(A, :, i))^2 for i in 1:size(A, 2)) + +function Base.LinAlg.qrfact(A::SparseMatrixCSC{Tv}; tol = _default_tol(A)) where {Tv <: CHOLMOD.VTypes} + R = Ref{Ptr{CHOLMOD.C_Sparse{Tv}}}() + E = Ref{Ptr{CHOLMOD.SuiteSparse_long}}() + H = Ref{Ptr{CHOLMOD.C_Sparse{Tv}}}() + HPinv = Ref{Ptr{CHOLMOD.SuiteSparse_long}}() + HTau = Ref{Ptr{CHOLMOD.C_Dense{Tv}}}(C_NULL) + + # SPQR doesn't accept symmetric matrices so we explicitly set the stype + r, p, hpinv = _qr!(ORDERING_DEFAULT, tol, 0, 0, Sparse(A, 0), + C_NULL, C_NULL, C_NULL, C_NULL, + R, E, H, HPinv, HTau) + + # convert to Julia managed memory and free memory from SuiteSparse + tmpH = Sparse(H[]) + HCSC = SparseMatrixCSC(tmpH) + free!(tmpH) + + tmpHTau = CHOLMOD.Dense(HTau[]) + HTauVector = vec(Array(tmpHTau)) + free!(tmpHTau) + + tmpR = Sparse(R[]) + RCSC = SparseMatrixCSC(tmpR) + free!(tmpR) + + return QRSparse(HCSC, HTauVector, RCSC, p, hpinv) end -function factorize(ordering::Integer, tol::Real, A::Sparse{Tv}) where Tv<:VTypes - s = unsafe_load(A.p) - if s.stype != 0 - throw(ArgumentError("stype must be zero")) +""" + qrfact(A) -> QRSparse + +Compute the `QR` factorization of a sparse matrix `A`. Fill-reducing row and column permutations +are used such that `F[:R] = F[:Q]'*A[F[:prow],F[:pcol]]`. The main application of this type is to +solve least squares or underdetermined problems with [`\\`](@ref). The function calls the C library SPQR. + +# Examples +```jldoctest +julia> A = sparse([1,2,3,4], [1,1,2,2], ones(4)) +4×2 SparseMatrixCSC{Float64,Int64} with 4 stored entries: + [1, 1] = 1.0 + [2, 1] = 1.0 + [3, 2] = 1.0 + [4, 2] = 1.0 + +julia> qrfact(A) +Base.SparseArrays.SPQR.QRSparse{Float64,Int64} +Q factor: +4×4 Base.SparseArrays.SPQR.QRSparseQ{Float64,Int64}: + -0.707107 0.0 0.0 -0.707107 + 0.0 -0.707107 -0.707107 0.0 + 0.0 -0.707107 0.707107 0.0 + -0.707107 0.0 0.0 0.707107 +R factor: +2×2 SparseMatrixCSC{Float64,Int64} with 2 stored entries: + [1, 1] = -1.41421 + [2, 2] = -1.41421 +Row permutation: +4-element Array{Int64,1}: + 1 + 3 + 4 + 2 +Columns permutation: +2-element Array{Int64,1}: + 1 + 2 +``` +""" +Base.LinAlg.qrfact(A::SparseMatrixCSC; tol = _default_tol(A)) = qrfact(A, Val{true}, tol = tol) + +Base.LinAlg.qr(A::SparseMatrixCSC; tol = _default_tol(A)) = qr(A, Val{true}, tol = tol) + +function Base.A_mul_B!(Q::QRSparseQ, A::StridedVecOrMat) + if size(A, 1) != size(Q, 1) + throw(DimensionMismatch("size(Q) = $(size(Q)) but size(A) = $(size(A))")) end - f = Factorization(size(A)..., ccall((:SuiteSparseQR_C_factorize, :libspqr), Ptr{C_Factorization{Tv}}, - (Cint, Cdouble, Ptr{Sparse{Tv}}, Ptr{Void}), - ordering, tol, get(A.p), common())) - finalizer(f, free!) - f + for l in size(Q.factors, 2):-1:1 + τl = -Q.τ[l] + h = view(Q.factors, :, l) + for j in 1:size(A, 2) + a = view(A, :, j) + LinAlg.axpy!(τl*dot(h, a), h, a) + end + end + return A end -function solve(system::Integer, QR::Factorization{Tv}, B::Dense{Tv}) where Tv<:VTypes - m, n = size(QR) - mB = size(B, 1) - if (system == RX_EQUALS_B || system == RETX_EQUALS_B) && m != mB - throw(DimensionMismatch("number of rows in factorized matrix must equal number of rows in right hand side")) - elseif (system == RTX_EQUALS_ETB || system == RTX_EQUALS_B) && n != mB - throw(DimensionMismatch("number of columns in factorized matrix must equal number of rows in right hand side")) +function Base.A_mul_B!(A::StridedMatrix, Q::QRSparseQ) + if size(A, 2) != size(Q, 1) + throw(DimensionMismatch("size(Q) = $(size(Q)) but size(A) = $(size(A))")) + end + tmp = similar(A, size(A, 1)) + for l in 1:size(Q.factors, 2) + τl = -Q.τ[l] + h = view(Q.factors, :, l) + A_mul_B!(tmp, A, h) + LinAlg.lowrankupdate!(A, tmp, h, τl) end - d = Dense(ccall((:SuiteSparseQR_C_solve, :libspqr), Ptr{C_Dense{Tv}}, - (Cint, Ptr{C_Factorization{Tv}}, Ptr{C_Dense{Tv}}, Ptr{Void}), - system, get(QR.p), get(B.p), common())) - finalizer(d, free!) - d + return A end -function qmult(method::Integer, QR::Factorization{Tv}, X::Dense{Tv}) where Tv<:VTypes - mQR, nQR = size(QR) - mX, nX = size(X) - if (method == QTX || method == QX) && mQR != mX - throw(DimensionMismatch("Q matrix size $mQR and dense matrix has $mX rows")) - elseif (method == XQT || method == XQ) && mQR != nX - throw(DimensionMismatch("Q matrix size $mQR and dense matrix has $nX columns")) +function Base.Ac_mul_B!(Q::QRSparseQ, A::StridedVecOrMat) + if size(A, 1) != size(Q, 1) + throw(DimensionMismatch("size(Q) = $(size(Q)) but size(A) = $(size(A))")) end - d = Dense(ccall((:SuiteSparseQR_C_qmult, :libspqr), Ptr{C_Dense{Tv}}, - (Cint, Ptr{C_Factorization{Tv}}, Ptr{C_Dense{Tv}}, Ptr{Void}), - method, get(QR.p), get(X.p), common())) - finalizer(d, free!) - d + for l in 1:size(Q.factors, 2) + τl = -Q.τ[l] + h = view(Q.factors, :, l) + for j in 1:size(A, 2) + a = view(A, :, j) + LinAlg.axpy!(τl'*dot(h, a), h, a) + end + end + return A end +function Base.A_mul_Bc!(A::StridedMatrix, Q::QRSparseQ) + if size(A, 2) != size(Q, 1) + throw(DimensionMismatch("size(Q) = $(size(Q)) but size(A) = $(size(A))")) + end + tmp = similar(A, size(A, 1)) + for l in size(Q.factors, 2):-1:1 + τl = -Q.τ[l] + h = view(Q.factors, :, l) + A_mul_B!(tmp, A, h) + LinAlg.lowrankupdate!(A, tmp, h, τl') + end + return A +end -qrfact(A::SparseMatrixCSC, ::Val{true}) = factorize(ORDERING_DEFAULT, DEFAULT_TOL, Sparse(A, 0)) +Base.LinAlg.getq(F::QRSparse) = QRSparseQ(F.factors, F.τ) """ - qrfact(A) -> SPQR.Factorization + getindex(F::QRSparse, d::Symbol) + +Extract factors of a QRSparse factorization. Possible values of `d` are +- `:Q` : `QRSparseQ` matrix of the ``Q`` factor in Householder form +- `:R` : `UpperTriangular` ``R`` factor +- `:prow` : Vector of the row permutations applied to the factorized matrix +- `:pcol` : Vector of the column permutations applied to the factorized matrix + +# Examples +```jldoctest +julia> F = qrfact(sparse([1,3,2,3,4], [1,1,2,3,4], [1.0,2.0,3.0,4.0,5.0])); + +julia> F[:Q] +4×4 Base.SparseArrays.SPQR.QRSparseQ{Float64,Int64}: + 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 + +julia> F[:R] +4×4 SparseMatrixCSC{Float64,Int64} with 5 stored entries: + [1, 1] = 3.0 + [2, 2] = 4.0 + [3, 3] = 5.0 + [2, 4] = 2.0 + [4, 4] = 1.0 -Compute the `QR` factorization of a sparse matrix `A`. A fill-reducing permutation is used. -The main application of this type is to solve least squares problems with [`\\`](@ref). The function -calls the C library SPQR and a few additional functions from the library are wrapped but not -exported. +julia> F[:prow] +4-element Array{Int64,1}: + 2 + 3 + 4 + 1 + +julia> F[:pcol] +4-element Array{Int64,1}: + 2 + 3 + 4 + 1 +``` """ -qrfact(A::SparseMatrixCSC) = qrfact(A, Val(true)) +function Base.getindex(F::QRSparse, d::Symbol) + if d == :Q + return LinAlg.getq(F) + elseif d == :R + return F.R + elseif d == :prow + return invperm(F.rpivinv) + elseif d == :pcol + return F.cpiv + else + throw(KeyError(d)) + end +end + +function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, F::QRSparse) + println(io, summary(F)) + println(io, "Q factor:") + show(io, mime, F[:Q]) + println(io, "\nR factor:") + show(io, mime, F[:R]) + println(io, "\nRow permutation:") + show(io, mime, F[:prow]) + println(io, "\nColumn permutation:") + show(io, mime, F[:pcol]) +end # 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 @@ -156,10 +343,10 @@ qrfact(A::SparseMatrixCSC) = qrfact(A, Val(true)) # here we have to use \ instead of A_ldiv_B! because of limitations in SPQR ## Two helper methods -_ret_size(F::Factorization, b::AbstractVector) = (size(F, 2),) -_ret_size(F::Factorization, B::AbstractMatrix) = (size(F, 2), size(B, 2)) +_ret_size(F::QRSparse, b::AbstractVector) = (size(F, 2),) +_ret_size(F::QRSparse, B::AbstractMatrix) = (size(F, 2), size(B, 2)) -function (\)(F::Factorization{Float64}, B::VecOrMat{Complex{Float64}}) +function (\)(F::QRSparse{Float64}, B::VecOrMat{Complex{Float64}}) # |z1|z3| reinterpret |x1|x2|x3|x4| transpose |x1|y1| reshape |x1|y1|x3|y3| # |z2|z4| -> |y1|y2|y3|y4| -> |x2|y2| -> |x2|y2|x4|y4| # |x3|y3| @@ -174,11 +361,65 @@ function (\)(F::Factorization{Float64}, B::VecOrMat{Complex{Float64}}) return reinterpret(Complex{Float64}, transpose(reshape(x, (length(x) >> 1), 2)), _ret_size(F, B)) end -function (\)(F::Factorization{T}, B::StridedVecOrMat{T}) where T<:VTypes - QtB = qmult(QTX, F, Dense(B)) - convert(typeof(B), solve(RETX_EQUALS_B, F, QtB)) +function _ldiv_basic(F::QRSparse, B::StridedVecOrMat) + if size(F, 1) != size(B, 1) + throw(DimensionMismatch("size(F) = $(size(F)) but size(B) = $(size(B))")) + end + + # The rank of F equal to the number of rows in R + rnk = size(F.R, 1) + + # allocate an array for the return value large enough to hold B and X + # For overdetermined problem, B is larger than X and vice versa + X = similar(B, ntuple(i -> i == 1 ? max(size(F, 2), size(B, 1)) : size(B, 2), Val(ndims(B)))) + + # Fill will zeros. These will eventually become the zeros in the basic solution + # fill!(X, 0) + # Apply left permutation to the solution and store in X + for j in 1:size(B, 2) + for i in 1:length(F.rpivinv) + @inbounds X[F.rpivinv[i], j] = B[i, j] + end + end + + # Make a view into X corresponding to the size of B + X0 = view(X, 1:size(B, 1), :) + + # Apply Q' to B + Ac_mul_B!(LinAlg.getq(F), X0) + + # Zero out to get basic solution + X[rnk + 1:end, :] = 0 + + # Solve R*X = B + A_ldiv_B!(UpperTriangular(view(F.R, :, Base.OneTo(rnk))), view(X0, Base.OneTo(rnk), :)) + + # Apply right permutation and extract solution from X + return getindex(X, ntuple(i -> i == 1 ? invperm(F.cpiv) : :, Val(ndims(B)))...) end -(\)(F::Factorization, B::StridedVecOrMat) = F\convert(AbstractArray{eltype(F)}, B) +(\)(F::QRSparse{T}, B::StridedVecOrMat{T}) where {T} = _ldiv_basic(F, B) +""" + (\)(F::QRSparse, B::StridedVecOrMat) + +Solve the least squares problem ``\min\|Ax - b\|^2`` or the linear system of equations +``Ax=b`` when `F` is the sparse QR factorization of ``A``. A basic solution is returned +when the problem is underdetermined. + +# Examples +```jldoctest +julia> A = sparse([1,2,4], [1,1,1], ones(3), 4, 2) +4×2 SparseMatrixCSC{Float64,Int64} with 3 stored entries: + [1, 1] = 1.0 + [2, 1] = 1.0 + [4, 1] = 1.0 + +julia> qrfact(A)\ones(4) +2-element Array{Float64,1}: + 1.0 + 0.0 +``` +""" +(\)(F::QRSparse, B::StridedVecOrMat) = F\convert(AbstractArray{eltype(F)}, B) end # module diff --git a/test/sparse/spqr.jl b/test/sparse/spqr.jl index eaa659dadd65c..0e17f0be2a065 100644 --- a/test/sparse/spqr.jl +++ b/test/sparse/spqr.jl @@ -3,17 +3,50 @@ using Base.SparseArrays.SPQR using Base.SparseArrays.CHOLMOD -let +@testset "Sparse QR" begin m, n = 100, 10 nn = 100 -for eltyA in (Float64, Complex{Float64}) - for eltyB in (Int, Float64, Complex{Float64}) - if eltyA <: Real - A = sparse([1:n; rand(1:m, nn - n)], [1:n; rand(1:n, nn - n)], randn(nn), m, n) - else - A = sparse([1:n; rand(1:m, nn - n)], [1:n; rand(1:n, nn - n)], complex.(randn(nn), randn(nn)), m, n) - end +@test size(qrfact(sprandn(m, n, 0.1))[:Q]) == (m, m) + +@testset "element type of A: $eltyA" for eltyA in (Float64, Complex{Float64}) + if eltyA <: Real + A = sparse([1:n; rand(1:m, nn - n)], [1:n; rand(1:n, nn - n)], randn(nn), m, n) + else + A = sparse([1:n; rand(1:m, nn - n)], [1:n; rand(1:n, nn - n)], complex.(randn(nn), randn(nn)), m, n) + end + + F = qrfact(A) + @test size(F) == (m,n) + @test size(F, 1) == m + @test size(F, 2) == n + @test size(F, 3) == 1 + @test_throws ArgumentError size(F, 0) + + @testset "getindex" begin + @test istriu(F[:R]) + @test isperm(F[:pcol]) + @test isperm(F[:prow]) + @test_throws KeyError F[:T] + end + + @testset "apply Q" begin + Q = F[:Q] + @test Q'*(Q*eye(m)) ≈ eye(m) + @test (eye(m)*Q)*Q' ≈ eye(m) + + # test that Q'Pl*A*Pr = R + R0 = Q'*full(A[F[:prow], F[:pcol]]) + @test R0[1:n, :] ≈ F[:R] + @test norm(R0[n + 1:end, :], 1) < 1e-12 + + @test_throws DimensionMismatch A_mul_B!(Q, eye(m + 1)) + @test_throws DimensionMismatch Ac_mul_B!(Q, eye(m + 1)) + @test_throws DimensionMismatch A_mul_B!(eye(m + 1), Q) + @test_throws DimensionMismatch A_mul_Bc!(eye(m + 1), Q) + end + + @testset "element type of B: $eltyB" for eltyB in (Int, Float64, Complex{Float64}) if eltyB == Int B = rand(1:10, m, 2) elseif eltyB <: Real @@ -27,34 +60,21 @@ for eltyA in (Float64, Complex{Float64}) @test A\B ≈ Array(A)\B @test_throws DimensionMismatch A\B[1:m-1,:] @test A[1:9,:]*(A[1:9,:]\ones(eltyB, 9)) ≈ ones(9) # Underdetermined system - - if eltyA == eltyB # promotions not defined for unexported methods - @test qrfact(sparse(eye(eltyA, 5)))\ones(eltyA, 5) == ones(5) - @test_throws ArgumentError SPQR.factorize(SPQR.ORDERING_DEFAULT, SPQR.DEFAULT_TOL, CHOLMOD.Sparse(sparse(eye(eltyA, 5)))) - @test_throws ArgumentError SPQR.Factorization(1, 1, convert(Ptr{SPQR.C_Factorization{eltyA}}, C_NULL)) - F = qrfact(A) - @test size(F) == (m,n) - @test size(F, 1) == m - @test size(F, 2) == n - @test size(F, 3) == 1 - @test_throws ArgumentError size(F, 0) - - # low level wrappers - @test_throws DimensionMismatch SPQR.solve(SPQR.RX_EQUALS_B, F, CHOLMOD.Dense(B')) - @test_throws DimensionMismatch SPQR.solve(SPQR.RTX_EQUALS_B, F, CHOLMOD.Dense(B)) - @test_throws DimensionMismatch SPQR.qmult(SPQR.QX, F, CHOLMOD.Dense(B')) - @test_throws DimensionMismatch SPQR.qmult(SPQR.XQ, F, CHOLMOD.Dense(B)) - @test A\B ≈ SPQR.backslash(SPQR.ORDERING_DEFAULT, SPQR.DEFAULT_TOL, CHOLMOD.Sparse(A), CHOLMOD.Dense(B)) - @test_throws DimensionMismatch SPQR.backslash(SPQR.ORDERING_DEFAULT, SPQR.DEFAULT_TOL, CHOLMOD.Sparse(A), CHOLMOD.Dense(B[1:m-1,:])) - end end + + # Make sure that conversion to Sparse doesn't use SuiteSparse's symmetric flag + @test qrfact(sparse(eye(eltyA, 5)))\ones(eltyA, 5) == ones(5) end -# Issue 14134 -F = qrfact(sprandn(10,5,0.5)) -b = IOBuffer() -serialize(b, F) -seekstart(b) -@test_throws ArgumentError deserialize(b)\ones(10) +@testset "basic solution of rank deficient ls" begin + A = sprandn(m, 5, 0.9)*sprandn(5, n, 0.9) + b = randn(m) + xs = A\b + xd = full(A)\b + + # check that basic solution has more zeros + @test countnz(xs) < countnz(xd) + @test A*xs ≈ A*xd +end end