Add qrfact(SparseMatrixCSC) by wrapping SPQR. Make \(SparseMatrixCSC)…

… work for least squares problems. Some tests.

base/linalg/generic.jl

@@ -237,11 +237,15 @@ function inv{T}(A::AbstractMatrix{T})
     A_ldiv_B!(factorize(convert(AbstractMatrix{S}, A)), eye(S, chksquare(A)))
 end
 
+function \{T}(A::AbstractMatrix{T}, B::AbstractVecOrMat{T})
+    size(A,1) == size(B,1) || throw(DimensionMismatch("LHS and RHS should have the same number of rows. LHS has $(size(A,1)) rows, but RHS has $(size(B,1)) rows."))
+    factorize(A)\B
+end
 function \{TA,TB}(A::AbstractMatrix{TA}, B::AbstractVecOrMat{TB})
     TC = typeof(one(TA)/one(TB))
-    size(A,1) == size(B,1) || throw(DimensionMismatch("LHS and RHS should have the same number of rows. LHS has $(size(A,1)) rows, but RHS has $(size(B,1)) rows."))
-    \(factorize(TA == TC ? A : convert(AbstractMatrix{TC}, A)), TB == TC ? copy(B) : convert(AbstractArray{TC}, B))
+    convert(AbstractMatrix{TC}, A)\convert(AbstractArray{TC}, B)
 end
+
 \(a::AbstractVector, b::AbstractArray) = reshape(a, length(a), 1) \ b
 /(A::AbstractVecOrMat, B::AbstractVecOrMat) = (B' \ A')'
 # \(A::StridedMatrix,x::Number) = inv(A)*x Should be added at some point when the old elementwise version has been deprecated long enough

base/sparse.jl

@@ -19,5 +19,6 @@ include("sparse/csparse.jl")
 include("sparse/linalg.jl")
 include("sparse/umfpack.jl")
 include("sparse/cholmod.jl")
+include("sparse/spqr.jl")
 
 end # module SparseMatrix

base/sparse/cholmod.jl

@@ -13,7 +13,7 @@ export
     Factor,
     Sparse
 
-using Base.SparseMatrix: AbstractSparseMatrix, SparseMatrixCSC, increment, increment!, indtype, decrement, decrement!
+using Base.SparseMatrix: AbstractSparseMatrix, SparseMatrixCSC, increment, indtype
 
 #########
 # Setup #
@@ -716,7 +716,7 @@ function Sparse(filename::ASCIIString)
 end
 
 ## convertion back to base Julia types
-function convert{T<:VTypes}(::Type{Matrix}, D::Dense{T})
+function convert{T}(::Type{Matrix{T}}, D::Dense{T})
     s = unsafe_load(D.p)
     a = Array(T, s.nrow, s.ncol)
     if s.d == s.nrow
@@ -730,6 +730,14 @@ function convert{T<:VTypes}(::Type{Matrix}, D::Dense{T})
     end
     a
 end
+convert{T}(::Type{Matrix}, D::Dense{T}) = convert(Matrix{T}, D)
+function convert{T}(::Type{Vector{T}}, D::Dense{T})
+    if size(D, 2) > 1
+        throw(DimensionMismatch("input must be a vector but had $(size(D, 2)) columns"))
+    end
+    reshape(convert(Matrix, D), size(D, 1))
+end
+convert{T}(::Type{Vector}, D::Dense{T}) = convert(Vector{T}, D)
 
 function convert{Tv,Ti}(::Type{SparseMatrixCSC{Tv,Ti}}, A::Sparse{Tv,Ti})
     s = unsafe_load(A.p)

base/sparse/linalg.jl

@@ -686,8 +686,10 @@ function factorize(A::SparseMatrixCSC)
             end
         end
         return lufact(A)
+    elseif m > n
+        return qrfact(A)
     else
-        throw(ArgumentError("sparse least squares problems by QR are not handled yet"))
+        throw(ArgumentError("underdetermined systemes are not implemented yet"))
     end
 end
 

base/sparse/spqr.jl

@@ -0,0 +1,133 @@
+module SPQR
+
+# ordering options */
+const ORDERING_FIXED   = int32(0)
+const ORDERING_NATURAL = int32(1)
+const ORDERING_COLAMD  = int32(2)
+const ORDERING_GIVEN   = int32(3) # only used for C/C++ interface
+const ORDERING_CHOLMOD = int32(4) # CHOLMOD best-effort (COLAMD, METIS,...)
+const ORDERING_AMD     = int32(5) # AMD(A'*A)
+const ORDERING_METIS   = int32(6) # metis(A'*A)
+const ORDERING_DEFAULT = int32(7) # SuiteSparseQR default ordering
+const ORDERING_BEST    = int32(8) # try COLAMD, AMD, and METIS; pick best
+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.
+
+# 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)
+
+# 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)
+
+
+using Base.SparseMatrix: SparseMatrixCSC
+using Base.SparseMatrix.CHOLMOD: C_Dense, C_Sparse, Dense, ITypes, Sparse, VTypes, common
+
+import Base: size
+import Base.LinAlg: qrfact
+import Base.SparseMatrix.CHOLMOD: convert, free!
+
+
+
+immutable C_Factorization{Tv<:VTypes,Ti<:ITypes}
+    xtype::Cint
+    factors::Ptr{Tv}
+end
+
+type Factorization{Tv<:VTypes,Ti<:ITypes} <: Base.LinAlg.Factorization{Tv}
+    m::Int
+    n::Int
+    p::Ptr{C_Factorization{Tv,Ti}}
+end
+
+size(F::Factorization) = (F.m, F.n)
+function size(F::Factorization, i::Integer)
+    if i < 1
+        throw(ArgumentError("dimension must be positive"))
+    end
+    if i == 1
+        return F.m
+    elseif i == 2
+        return F.n
+    end
+    return 1
+end
+
+function free!{Tv<:VTypes,Ti<:ITypes}(F::Factorization{Tv,Ti})
+    ccall((:SuiteSparseQR_C_free, :libspqr), Cint,
+        (Ptr{Ptr{C_Factorization{Tv,Ti}}}, Ptr{Void}),
+            &F.p, common(Ti)) == 1
+end
+
+function backslash{Tv<:VTypes,Ti<:ITypes}(ordering::Integer, tol::Real, A::Sparse{Tv,Ti}, B::Dense{Tv})
+    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,Ti}}, Ptr{C_Dense{Tv}}, Ptr{Void}),
+            ordering, tol, A.p, B.p, common(Ti)))
+    finalizer(d, free!)
+    d
+end
+
+function factorize{Tv<:VTypes,Ti<:ITypes}(ordering::Integer, tol::Real, A::Sparse{Tv,Ti})
+    f = Factorization(size(A)..., ccall((:SuiteSparseQR_C_factorize, :libspqr), Ptr{C_Factorization{Tv,Ti}},
+        (Cint, Cdouble, Ptr{Sparse{Tv,Ti}}, Ptr{Void}),
+            ordering, tol, A.p, common(Ti)))
+    finalizer(f, free!)
+    f
+end
+
+function solve{Tv<:VTypes,Ti<:ITypes}(system::Integer, QR::Factorization{Tv,Ti}, B::Dense{Tv})
+    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"))
+    end
+    d = Dense(ccall((:SuiteSparseQR_C_solve, :libspqr), Ptr{C_Dense{Tv}},
+        (Cint, Ptr{C_Factorization{Tv,Ti}}, Ptr{C_Dense{Tv}}, Ptr{Void}),
+            system, QR.p, B.p, common(Ti)))
+    finalizer(d, free!)
+    d
+end
+
+function qmult{Tv<:VTypes,Ti<:ITypes}(method::Integer, QR::Factorization{Tv,Ti}, X::Dense{Tv})
+    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"))
+    end
+    d = Dense(ccall((:SuiteSparseQR_C_qmult, :libspqr), Ptr{C_Dense{Tv}},
+            (Cint, Ptr{C_Factorization{Tv,Ti}}, Ptr{C_Dense{Tv}}, Ptr{Void}),
+                method, QR.p, X.p, common(Ti)))
+    finalizer(d, free!)
+    d
+end
+
+qrfact(A::SparseMatrixCSC) = factorize(ORDERING_DEFAULT, DEFAULT_TOL, Sparse(A))
+
+function (\){T}(F::Factorization{T}, B::StridedVecOrMat{T})
+    QtB = qmult(QTX, F, Dense(B))
+    convert(typeof(B), solve(RETX_EQUALS_B, F, QtB))
+end
+
+end # module

test/runtests.jl

@@ -3,7 +3,7 @@ testnames = [
     "linalg", "core", "keywordargs", "numbers", "strings", "dates",
     "dict", "hashing", "remote", "iobuffer", "staged", "arrayops",
     "subarray", "reduce", "reducedim", "random", "intfuncs",
-    "simdloop", "blas", "fft", "dsp", "sparse", "bitarray", "copy", "math",
+    "simdloop", "blas", "fft", "dsp", "sparsetests", "bitarray", "copy", "math",
     "fastmath", "functional", "bigint", "sorting", "statistics", "spawn",
     "backtrace", "priorityqueue", "arpack", "file", "version",
     "resolve", "pollfd", "mpfr", "broadcast", "complex", "socket",
@@ -25,12 +25,6 @@ push!(testnames, "parallel")
 
 tests = (ARGS==["all"] || isempty(ARGS)) ? testnames : ARGS
 
-if "sparse" in tests
-    # specifically selected case
-    filter!(x -> x != "sparse", tests)
-    prepend!(tests, ["sparse/sparse", "sparse/cholmod", "sparse/umfpack"])
-end
-
 if "linalg" in tests
     # specifically selected case
     filter!(x -> x != "linalg", tests)

test/sparse/spqr.jl

@@ -0,0 +1,43 @@
+using Base.Test
+
+using Base.SparseMatrix.SPQR
+using Base.SparseMatrix.CHOLMOD
+
+m, n = 1000, 10
+nn = 100
+
+for eltyA in (Float64, Complex{Float64})
+    for eltyB in (Float64, Complex{Float64})
+        if eltyA <: Real
+            A = sparse(rand(1:m, nn), rand(1:n, nn), randn(nn), m, n)
+        else
+            A = sparse(rand(1:m, nn), rand(1:n, nn), complex(randn(nn), randn(nn)), m, n)
+        end
+        if eltyB <: Real
+            B = randn(m, 2)
+        else
+            B = complex(randn(m, 2), randn(m, 2))
+        end
+
+        @test_approx_eq A\B[:,1] full(A)\B[:,1]
+        @test_approx_eq A\B full(A)\B
+        @test_throws DimensionMismatch A\B[1:m-1,:]
+
+        if eltyA == eltyB # promotions not defined for unexported methods
+            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_approx_eq 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
+end

test/sparsetests.jl

@@ -0,0 +1,4 @@
+include("sparse/sparse.jl")
+include("sparse/umfpack.jl")
+include("sparse/cholmod.jl")
+include("sparse/spqr.jl")