Skip to content

A[c|t]_ldiv_B! specializations for UmfpackLU-StridedVecOrMat, less generalized linear indexing and meta-fu, take 2 #20055

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 1 commit into from
Jan 27, 2017
Merged

A[c|t]_ldiv_B! specializations for UmfpackLU-StridedVecOrMat, less generalized linear indexing and meta-fu, take 2 #20055

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
87 changes: 53 additions & 34 deletions base/sparse/umfpack.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -383,42 +383,61 @@ function nnz(lu::UmfpackLU)
end

### Solve with Factorization
for (f!, umfpack) in ((:A_ldiv_B!, :UMFPACK_A),
(:Ac_ldiv_B!, :UMFPACK_At),
(:At_ldiv_B!, :UMFPACK_Aat))
@eval begin
function $f!{T<:UMFVTypes}(x::StridedVecOrMat{T}, lu::UmfpackLU{T}, b::StridedVecOrMat{T})
n = size(x, 2)
if n != size(b, 2)
throw(DimensionMismatch("in and output arrays must have the same number of columns"))
end
for j in 1:n
solve!(view(x, :, j), lu, view(b, :, j), $umfpack)
end
return x
end
$f!{T<:UMFVTypes}(lu::UmfpackLU{T}, b::StridedVector{T}) = $f!(b, lu, copy(b))
$f!{T<:UMFVTypes}(lu::UmfpackLU{T}, b::StridedMatrix{T}) = $f!(b, lu, copy(b))

function $f!{Tb<:Complex}(x::StridedVector{Tb}, lu::UmfpackLU{Float64}, b::StridedVector{Tb})
m, n = size(x, 1), size(x, 2)
if n != size(b, 2)
throw(DimensionMismatch("in and output arrays must have the same number of columns"))
end
# TODO: Optionally let user allocate these and pass in somehow
r = similar(b, Float64, m)
i = similar(b, Float64, m)
for j in 1:n
solve!(r, lu, convert(Vector{Float64}, real(view(b, :, j))), $umfpack)
solve!(i, lu, convert(Vector{Float64}, imag(view(b, :, j))), $umfpack)

map!((t,s) -> t + im*s, view(x, :, j), r, i)
end
return x
end
$f!{Tb<:Complex}(lu::UmfpackLU{Float64}, b::StridedVector{Tb}) = $f!(b, lu, copy(b))
A_ldiv_B!{T<:UMFVTypes}(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) = A_ldiv_B!(B, lu, copy(B))
At_ldiv_B!{T<:UMFVTypes}(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) = At_ldiv_B!(B, lu, copy(B))
Ac_ldiv_B!{T<:UMFVTypes}(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) = Ac_ldiv_B!(B, lu, copy(B))
A_ldiv_B!{Tb<:Complex}(lu::UmfpackLU{Float64}, B::StridedVecOrMat{Tb}) = A_ldiv_B!(B, lu, copy(B))
At_ldiv_B!{Tb<:Complex}(lu::UmfpackLU{Float64}, B::StridedVecOrMat{Tb}) = At_ldiv_B!(B, lu, copy(B))
Ac_ldiv_B!{Tb<:Complex}(lu::UmfpackLU{Float64}, B::StridedVecOrMat{Tb}) = Ac_ldiv_B!(B, lu, copy(B))

A_ldiv_B!{T<:UMFVTypes}(X::StridedVecOrMat{T}, lu::UmfpackLU{T}, B::StridedVecOrMat{T}) =
_Aq_ldiv_B!(X, lu, B, UMFPACK_A)
At_ldiv_B!{T<:UMFVTypes}(X::StridedVecOrMat{T}, lu::UmfpackLU{T}, B::StridedVecOrMat{T}) =
_Aq_ldiv_B!(X, lu, B, UMFPACK_At)
Ac_ldiv_B!{T<:UMFVTypes}(X::StridedVecOrMat{T}, lu::UmfpackLU{T}, B::StridedVecOrMat{T}) =
_Aq_ldiv_B!(X, lu, B, UMFPACK_Aat)
A_ldiv_B!{Tb<:Complex}(X::StridedVecOrMat{Tb}, lu::UmfpackLU{Float64}, B::StridedVecOrMat{Tb}) =
_Aq_ldiv_B!(X, lu, B, UMFPACK_A)
At_ldiv_B!{Tb<:Complex}(X::StridedVecOrMat{Tb}, lu::UmfpackLU{Float64}, B::StridedVecOrMat{Tb}) =
_Aq_ldiv_B!(X, lu, B, UMFPACK_At)
Ac_ldiv_B!{Tb<:Complex}(X::StridedVecOrMat{Tb}, lu::UmfpackLU{Float64}, B::StridedVecOrMat{Tb}) =
Copy link
Contributor

@tkelman tkelman Jan 27, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What happens here if one input's a vector and the other a 1-column matrix? I think that would have worked before but doesn't now?

Copy link
Member Author

@Sacha0 Sacha0 Feb 18, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good catch! I believe you are correct. From brief inspection, covering the vector-colmat and colmat-vector cases would require four additional kernels. Would covering those cases justify the additional code complexity? Thanks!

Copy link
Contributor

@tkelman tkelman Feb 18, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

if that functionality was being used anywhere, then might be a reason to prefer the old version

Copy link
Member Author

@Sacha0 Sacha0 Feb 18, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Of course please feel welcome to revert this commit if so. (And if so, we should revive #20045 to retain the bug fix.) Thanks!

_Aq_ldiv_B!(X, lu, B, UMFPACK_Aat)

function _Aq_ldiv_B!(X::StridedVecOrMat, lu::UmfpackLU, B::StridedVecOrMat, transposeoptype)
if size(X, 2) != size(B, 2)
throw(DimensionMismatch("input and output arrays must have same number of columns"))
end
_AqldivB_kernel!(X, lu, B, transposeoptype)
return X
end
function _AqldivB_kernel!{T<:UMFVTypes}(x::StridedVector{T}, lu::UmfpackLU{T},
b::StridedVector{T}, transposeoptype)
solve!(x, lu, b, transposeoptype)
end
function _AqldivB_kernel!{T<:UMFVTypes}(X::StridedMatrix{T}, lu::UmfpackLU{T},
B::StridedMatrix{T}, transposeoptype)
for col in 1:size(X, 2)
solve!(view(X, :, col), lu, view(B, :, col), transposeoptype)
end
end
function _AqldivB_kernel!{Tb<:Complex}(x::StridedVector{Tb}, lu::UmfpackLU{Float64},
b::StridedVector{Tb}, transposeoptype)
r, i = similar(b, Float64), similar(b, Float64)
solve!(r, lu, Vector{Float64}(real(b)), transposeoptype)
solve!(i, lu, Vector{Float64}(imag(b)), transposeoptype)
map!(complex, x, r, i)
end
function _AqldivB_kernel!{Tb<:Complex}(X::StridedMatrix{Tb}, lu::UmfpackLU{Float64},
B::StridedMatrix{Tb}, transposeoptype)
r = similar(B, Float64, size(B, 1))
i = similar(B, Float64, size(B, 1))
for j in 1:size(B, 2)
solve!(r, lu, Vector{Float64}(real(view(B, :, j))), transposeoptype)
solve!(i, lu, Vector{Float64}(imag(view(B, :, j))), transposeoptype)
map!(complex, view(X, :, j), r, i)
end
end


function getindex(lu::UmfpackLU, d::Symbol)
L,U,p,q,Rs = umf_extract(lu)
Expand Down
12 changes: 12 additions & 0 deletions test/sparse/umfpack.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -143,3 +143,15 @@ let
F = lufact(A)
@test F[:p] == [3 ; 4 ; 2 ; 1]
end

# Test that A[c|t]_ldiv_B!{T<:Complex}(X::StridedMatrix{T}, lu::UmfpackLU{Float64},
# B::StridedMatrix{T}) works as expected.
let N = 10, p = 0.5
A = N*speye(N) + sprand(N, N, p)
X = zeros(Complex{Float64}, N, N)
B = complex.(rand(N, N), rand(N, N))
luA, lufA = lufact(A), lufact(Array(A))
@test A_ldiv_B!(copy(X), luA, B) ≈ A_ldiv_B!(copy(X), lufA, B)
@test At_ldiv_B!(copy(X), luA, B) ≈ At_ldiv_B!(copy(X), lufA, B)
@test Ac_ldiv_B!(copy(X), luA, B) ≈ Ac_ldiv_B!(copy(X), lufA, B)
end