Merge pull request #13793 from JuliaLang/anj/sparsesolve

RFC: Use sparse triangular solvers for sparse triangular solves.

base/linalg/dense.jl

@@ -436,20 +436,6 @@ function factorize{T}(A::Matrix{T})
     qrfact(A, Val{true})
 end
 
-(\)(a::Vector, B::StridedVecOrMat) = (\)(reshape(a, length(a), 1), B)
-
-function (\)(A::StridedMatrix, B::StridedVecOrMat)
-    m, n = size(A)
-    if m == n
-        if istril(A)
-            return istriu(A) ? \(Diagonal(A),B) : \(LowerTriangular(A),B)
-        end
-        istriu(A) && return \(UpperTriangular(A),B)
-        return \(lufact(A),B)
-    end
-    return qrfact(A,Val{true})\B
-end
-
 ## Moore-Penrose pseudoinverse
 function pinv{T}(A::StridedMatrix{T}, tol::Real)
     m, n = size(A)

base/linalg/diagonal.jl

@@ -5,7 +5,7 @@
 immutable Diagonal{T} <: AbstractMatrix{T}
     diag::Vector{T}
 end
-Diagonal(A::Matrix) = Diagonal(diag(A))
+Diagonal(A::AbstractMatrix) = Diagonal(diag(A))
 
 convert{T}(::Type{Diagonal{T}}, D::Diagonal{T}) = D
 convert{T}(::Type{Diagonal{T}}, D::Diagonal) = Diagonal{T}(convert(Vector{T}, D.diag))
@@ -163,9 +163,9 @@ function A_ldiv_B!(D::Diagonal, B::StridedVecOrMat)
     end
     return B
 end
-\(D::Diagonal, B::StridedMatrix) = scale(1 ./ D.diag, B)
-\(D::Diagonal, b::StridedVector) = reshape(scale(1 ./ D.diag, reshape(b, length(b), 1)), length(b))
-\(Da::Diagonal, Db::Diagonal) = Diagonal(Db.diag ./ Da.diag)
+(\)(D::Diagonal, B::AbstractMatrix) = scale(1 ./ D.diag, B)
+(\)(D::Diagonal, b::AbstractVector) = reshape(scale(1 ./ D.diag, reshape(b, length(b), 1)), length(b))
+(\)(Da::Diagonal, Db::Diagonal) = Diagonal(Db.diag ./ Da.diag)
 
 function inv{T}(D::Diagonal{T})
     Di = similar(D.diag)

base/linalg/factorization.jl

@@ -20,7 +20,21 @@ end
 ### General promotion rules
 convert{T}(::Type{Factorization{T}}, F::Factorization{T}) = F
 inv{T}(F::Factorization{T}) = A_ldiv_B!(F, eye(T, size(F,1)))
-function \{TF<:Number,TB<:Number,N}(F::Factorization{TF}, B::AbstractArray{TB,N})
+
+# 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 (\){T<:BlasReal}(F::Factorization{T}, B::AbstractVector{Complex{T}})
+    c2r = reshape(transpose(reinterpret(T, B, (2, length(B)))), size(B, 1), 2*size(B, 2))
+    x = A_ldiv_B!(F, c2r)
+    return reinterpret(Complex{T}, transpose(reshape(x, div(length(x), 2), 2)), (size(F,2),))
+end
+function (\){T<:BlasReal}(F::Factorization{T}, B::AbstractMatrix{Complex{T}})
+    c2r = reshape(transpose(reinterpret(T, B, (2, length(B)))), size(B, 1), 2*size(B, 2))
+    x = A_ldiv_B!(F, c2r)
+    return reinterpret(Complex{T}, transpose(reshape(x, div(length(x), 2), 2)), (size(F,2), size(B,2)))
+end
+
+function (\){TF<:Number,TB<:Number,N}(F::Factorization{TF}, B::AbstractArray{TB,N})
     TFB = typeof(one(TF)/one(TB))
     A_ldiv_B!(convert(Factorization{TFB}, F), copy_oftype(B, TFB))
 end

base/linalg/generic.jl

@@ -325,19 +325,26 @@ 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})
-    if size(A,1) != size(B,1)
-        throw(DimensionMismatch("left and right hand sides should have the same number of rows, left hand side has $(size(A,1)) rows, but right hand side has $(size(B,1)) rows."))
+function (\)(A::AbstractMatrix, B::AbstractVecOrMat)
+    m, n = size(A)
+    if m == n
+        if istril(A)
+            if istriu(A)
+                return Diagonal(A) \ B
+            else
+                return LowerTriangular(A) \ B
+            end
+        end
+        if istriu(A)
+            return UpperTriangular(A) \ B
+        end
+        return lufact(A) \ B
     end
-    factorize(A)\B
-end
-function \{TA,TB}(A::AbstractMatrix{TA}, B::AbstractVecOrMat{TB})
-    TC = typeof(one(TA)/one(TB))
-    convert(AbstractMatrix{TC}, A)\convert(AbstractArray{TC}, B)
+    return qrfact(A,Val{true}) \ B
 end
 
-\(a::AbstractVector, b::AbstractArray) = reshape(a, length(a), 1) \ b
-/(A::AbstractVecOrMat, B::AbstractVecOrMat) = (B' \ A')'
+(\)(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
 # /(x::Number,A::StridedMatrix) = x*inv(A)
 

base/sparse.jl

@@ -6,7 +6,7 @@ using Base: Func, AddFun, OrFun, ConjFun, IdFun
 using Base.Sort: Forward
 using Base.LinAlg: AbstractTriangular, PosDefException
 
-import Base: +, -, *, &, |, $, .+, .-, .*, ./, .\, .^, .<, .!=, ==
+import Base: +, -, *, \, &, |, $, .+, .-, .*, ./, .\, .^, .<, .!=, ==
 import Base: A_mul_B!, Ac_mul_B, Ac_mul_B!, At_mul_B, At_mul_B!, A_ldiv_B!
 
 import Base: @get!, acos, acosd, acot, acotd, acsch, asech, asin, asind, asinh,

base/sparse/linalg.jl

@@ -167,16 +167,10 @@ end
 ## solvers
 function fwdTriSolve!(A::SparseMatrixCSC, B::AbstractVecOrMat)
 # forward substitution for CSC matrices
-    n = length(B)
-    if isa(B, Vector)
-        nrowB = n
-        ncolB = 1
-    else
-        nrowB, ncolB = size(B)
-    end
-    ncol = chksquare(A)
+    nrowB, ncolB  = size(B, 1), size(B, 2)
+    ncol = LinAlg.chksquare(A)
     if nrowB != ncol
-        throw(DimensionMismatch("A is $(ncol)X$(ncol) and B has length $(n)"))
+        throw(DimensionMismatch("A is $(ncol) columns and B has $(nrowB) rows"))
     end
 
     aa = A.nzval
@@ -185,56 +179,88 @@ function fwdTriSolve!(A::SparseMatrixCSC, B::AbstractVecOrMat)
 
     joff = 0
     for k = 1:ncolB
-        for j = 1:(nrowB-1)
-            jb = joff + j
+        for j = 1:nrowB
             i1 = ia[j]
-            i2 = ia[j+1]-1
-            B[jb] /= aa[i1]
-            bj = B[jb]
-            for i = i1+1:i2
-                B[joff+ja[i]] -= bj*aa[i]
+            i2 = ia[j + 1] - 1
+
+            # loop through the structural zeros
+            ii = i1
+            jai = ja[ii]
+            while ii <= i2 && jai < j
+                ii += 1
+                jai = ja[ii]
+            end
+
+            # check for zero pivot and divide with pivot
+            if jai == j
+                bj = B[joff + jai]/aa[ii]
+                B[joff + jai] = bj
+                ii += 1
+            else
+                throw(LinAlg.SingularException(j))
+            end
+
+            # update remaining part
+            for i = ii:i2
+                B[joff + ja[i]] -= bj*aa[i]
             end
         end
         joff += nrowB
-        B[joff] /= aa[end]
     end
-    return B
+    B
 end
 
 function bwdTriSolve!(A::SparseMatrixCSC, B::AbstractVecOrMat)
 # backward substitution for CSC matrices
-    n = length(B)
-    if isa(B, Vector)
-        nrowB = n
-        ncolB = 1
-    else
-        nrowB, ncolB = size(B)
+    nrowB, ncolB = size(B, 1), size(B, 2)
+    ncol = LinAlg.chksquare(A)
+    if nrowB != ncol
+        throw(DimensionMismatch("A is $(ncol) columns and B has $(nrowB) rows"))
     end
-    ncol = chksquare(A)
-    if nrowB != ncol throw(DimensionMismatch("A is $(ncol)X$(ncol) and B has length $(n)")) end
 
     aa = A.nzval
     ja = A.rowval
     ia = A.colptr
 
     joff = 0
     for k = 1:ncolB
-        for j = nrowB:-1:2
-            jb = joff + j
+        for j = nrowB:-1:1
             i1 = ia[j]
-            i2 = ia[j+1]-1
-            B[jb] /= aa[i2]
-            bj = B[jb]
-            for i = i2-1:-1:i1
-                B[joff+ja[i]] -= bj*aa[i]
+            i2 = ia[j + 1] - 1
+
+            # loop through the structural zeros
+            ii = i2
+            jai = ja[ii]
+            while ii >= i1 && jai > j
+                ii -= 1
+                jai = ja[ii]
+            end
+
+            # check for zero pivot and divide with pivot
+            if jai == j
+                bj = B[joff + jai]/aa[ii]
+                B[joff + jai] = bj
+                ii -= 1
+            else
+                throw(LinAlg.SingularException(j))
+            end
+
+            # update remaining part
+            for i = ii:-1:i1
+                B[joff + ja[i]] -= bj*aa[i]
             end
         end
-        B[joff+1] /= aa[1]
         joff += nrowB
     end
-   return B
+    B
 end
 
+A_ldiv_B!{T,Ti}(L::LowerTriangular{T,SparseMatrixCSC{T,Ti}}, B::StridedVecOrMat) = fwdTriSolve!(L.data, B)
+A_ldiv_B!{T,Ti}(U::UpperTriangular{T,SparseMatrixCSC{T,Ti}}, B::StridedVecOrMat) = bwdTriSolve!(U.data, B)
+
+(\){T,Ti}(L::LowerTriangular{T,SparseMatrixCSC{T,Ti}}, B::SparseMatrixCSC) = A_ldiv_B!(L, full(B))
+(\){T,Ti}(U::UpperTriangular{T,SparseMatrixCSC{T,Ti}}, B::SparseMatrixCSC) = A_ldiv_B!(U, full(B))
+
 ## triu, tril
 
 function triu{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti}, k::Integer=0)
@@ -800,6 +826,15 @@ scale{T,Tv,Ti}(b::Vector{T}, A::SparseMatrixCSC{Tv,Ti}) =
 function factorize(A::SparseMatrixCSC)
     m, n = size(A)
     if m == n
+        if istril(A)
+            if istriu(A)
+                return return Diagonal(A)
+            else
+                return LowerTriangular(A)
+            end
+        elseif istriu(A)
+            return UpperTriangular(A)
+        end
         AC = CHOLMOD.Sparse(A)
         if ishermitian(AC)
             try

base/sparse/spqr.jl

@@ -138,12 +138,28 @@ function qmult{Tv<:VTypes}(method::Integer, QR::Factorization{Tv}, X::Dense{Tv})
     d
 end
 
-qrfact(A::SparseMatrixCSC) = factorize(ORDERING_DEFAULT, DEFAULT_TOL, Sparse(A, 0))
-
-function (\){T}(F::Factorization{T}, B::StridedVecOrMat{T})
+qrfact(A::SparseMatrixCSC, ::Type{Val{true}}) = factorize(ORDERING_DEFAULT, DEFAULT_TOL, Sparse(A, 0))
+qrfact(A::SparseMatrixCSC) = qrfact(A, Val{true})
+
+# 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
+#
+# This definition is similar to the definition in factorization.jl except the we here have to use
+# \ instead of A_ldiv_B! because of limitations in SPQR
+function (\)(F::Factorization{Float64}, B::StridedVector{Complex{Float64}})
+    c2r = reshape(transpose(reinterpret(Float64, B, (2, length(B)))), size(B, 1), 2*size(B, 2))
+    x = F\c2r
+    return reinterpret(Complex{Float64}, transpose(reshape(x, div(length(x), 2), 2)), (size(F,2),))
+end
+function (\)(F::Factorization{Float64}, B::StridedMatrix{Complex{Float64}})
+    c2r = reshape(transpose(reinterpret(Float64, B, (2, length(B)))), size(B, 1), 2*size(B, 2))
+    x = F\c2r
+    return reinterpret(Complex{Float64}, transpose(reshape(x, div(length(x), 2), 2)), (size(F,2), size(B,2)))
+end
+function (\){T<:VTypes}(F::Factorization{T}, B::StridedVecOrMat{T})
     QtB = qmult(QTX, F, Dense(B))
     convert(typeof(B), solve(RETX_EQUALS_B, F, QtB))
 end
-(\){T,S}(F::Factorization{T}, B::StridedVecOrMat{S}) = F\convert(AbstractArray{T}, B)
+(\)(F::Factorization, B::StridedVecOrMat) = F\convert(AbstractArray{eltype(F)}, B)
 
 end # module

base/sparse/umfpack.jl

@@ -115,9 +115,22 @@ function lufact{Tv<:UMFVTypes,Ti<:UMFITypes}(S::SparseMatrixCSC{Tv,Ti})
 end
 lufact(A::SparseMatrixCSC) = lufact(float(A))
 
-function show(io::IO, f::UmfpackLU)
-    println(io, "UMFPACK LU Factorization of a $(f.m)-by-$(f.n) sparse matrix")
-    f.numeric != C_NULL && println(io, f.numeric)
+size(F::UmfpackLU) = (F.m, F.n)
+function size(F::UmfpackLU, dim::Integer)
+    if dim < 1
+        error("arraysize: dimension out of range")
+    elseif dim == 1
+        return Int(F.m)
+    elseif dim == 2
+        return Int(F.n)
+    else
+        return 1
+    end
+end
+
+function show(io::IO, F::UmfpackLU)
+    println(io, "UMFPACK LU Factorization of a $(size(F)) sparse matrix")
+    F.numeric != C_NULL && println(io, F.numeric)
 end
 
 ## Wrappers for UMFPACK functions

test/sparsedir/cholmod.jl

@@ -174,7 +174,6 @@ let # Issue 9160
     end
 end
 
-
 # Issue #9915
 @test speye(2)\speye(2) == eye(2)
 

test/sparsedir/sparse.jl

@@ -1175,3 +1175,15 @@ let
     @test_throws ErrorException eig(A)
     @test_throws ErrorException inv(A)
 end
+
+let
+    n = 100
+    A = sprandn(n, n, 0.5) + sqrt(n)*I
+    x = LowerTriangular(A)*ones(n)
+    @test LowerTriangular(A)\x ≈ ones(n)
+    x = UpperTriangular(A)*ones(n)
+    @test UpperTriangular(A)\x ≈ ones(n)
+    A[2,2] = 0
+    @test_throws LinAlg.SingularException LowerTriangular(A)\ones(n)
+    @test_throws LinAlg.SingularException UpperTriangular(A)\ones(n)
+end