diff --git a/Project.toml b/Project.toml index 5281432..c2103e1 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "KLU" uuid = "ef3ab10e-7fda-4108-b977-705223b18434" authors = ["Raye Kimmerer and contributors"] -version = "0.6.1" +version = "0.7.0" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/README.md b/README.md index 9e5e972..1b172cc 100644 --- a/README.md +++ b/README.md @@ -37,6 +37,12 @@ julia> x = factor \ b 10.688731342850245 ``` +## Supported types + +KLU provides `Float64` and `ComplexF64` values with `Int32` and `Int64` indices, and nothing else is +accepted. Right-hand sides must have the same element type as the factorization. Convert anything +else first, for example with `klu(Float64.(A))` or `factor \ Float64.(b)`. + ## Reusing a factorization If only the values of a matrix change, `klu!` refactors faster than `klu`. It reuses the @@ -52,9 +58,9 @@ If `klu!` throws a `SingularException`, `factor` has no numeric factorization un `klu!` or `klu_factor!` succeeds, and solving with it throws as well. To avoid using more memory, `factor` shares the vectors it is given: after `klu(A)` or -`klu!(factor, B)` with `Float64` or `ComplexF64` values, `nonzeros(factor)` is the vector of `A` -or `B`. It is the caller's responsibility not to change them in a way that breaks `factor`, as -described in the docstring of `klu`. +`klu!(factor, B)`, `nonzeros(factor)` is the vector of `A` or `B`. It is the caller's +responsibility not to change them in a way that breaks `factor`, as described in the docstring of +`klu`. `klu(A; full_factor = false)` only analyzes the pattern of `A`, and `klu!` then computes the first numeric factorization. Use `copy(factor)` to get an independent factorization, for diff --git a/docs/src/index.md b/docs/src/index.md index f5c5d7f..b34becf 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -38,6 +38,12 @@ julia> x = factor \ b 10.688731342850245 ``` +## Supported types + +KLU provides `Float64` and `ComplexF64` values with `Int32` and `Int64` indices, and nothing else is +accepted. Right-hand sides must have the same element type as the factorization. Convert anything +else first, for example with `klu(Float64.(A))` or `factor \ Float64.(b)`. + ## Reusing a factorization If only the values of a matrix change, `klu!` refactors faster than `klu`. It reuses the @@ -53,9 +59,9 @@ If `klu!` throws a `SingularException`, `factor` has no numeric factorization un `klu!` or `klu_factor!` succeeds, and solving with it throws as well. To avoid using more memory, `factor` shares the vectors it is given: after `klu(A)` or -`klu!(factor, B)` with `Float64` or `ComplexF64` values, `nonzeros(factor)` is the vector of `A` -or `B`. It is the caller's responsibility not to change them in a way that breaks `factor`, as -described in the docstring of [`klu`](@ref). +`klu!(factor, B)`, `nonzeros(factor)` is the vector of `A` or `B`. It is the caller's +responsibility not to change them in a way that breaks `factor`, as described in the docstring of +[`klu`](@ref). `klu(A; full_factor = false)` only analyzes the pattern of `A`, and `klu!` then computes the first numeric factorization. Use `copy(factor)` to get an independent factorization, for diff --git a/src/KLU.jl b/src/KLU.jl index 5c3078a..4239ced 100644 --- a/src/KLU.jl +++ b/src/KLU.jl @@ -172,10 +172,9 @@ mutable struct KLUFactorization{Tv<:KLUTypes, Ti<:KLUITypes} <: AbstractKLUFacto const lock::ReentrantLock # The permutations given to `klu_analyze!`, if any, so that a copy can analyze with them. _perms::Union{Nothing, Tuple{Vector{Ti}, Vector{Ti}}} - function KLUFactorization(n, colptr, rowval, nzval, common = _common(eltype(colptr))) + function KLUFactorization(n, colptr::Vector{Ti}, rowval::Vector{Ti}, nzval::Vector{Tv}, common = _common(Ti)) where {Tv<:KLUTypes, Ti<:KLUITypes} _check_lengths(n, colptr, rowval, nzval) - Ti = eltype(colptr) - return new{eltype(nzval), Ti}(common, C_NULL, C_NULL, n, colptr, rowval, nzval, ReentrantLock(), nothing) + return new{Tv, Ti}(common, C_NULL, C_NULL, n, colptr, rowval, nzval, ReentrantLock(), nothing) end end @@ -246,9 +245,8 @@ nnz(K::AbstractKLUFactorization) = Base.@lock K K.lnz + K.unz + K.nzoff The values of the matrix factored by `K`. -!!! warning "nonzeros(K) may or may not point to nonzeros(A)" - After `K = klu(A)`, `nonzeros(K) === nonzeros(A)` if `A` has `Float64` or `ComplexF64` values, - so changing one changes the other. Other value types are converted to a new vector. +!!! warning "nonzeros(K) points to nonzeros(A)" + After `K = klu(A)`, `nonzeros(K) === nonzeros(A)`, so changing one changes the other. See [`klu`](@ref) for what changing these values means for `K`. """ nonzeros(K::AbstractKLUFactorization) = K.nzval @@ -703,10 +701,9 @@ klu_factor! Compute the LU factorization of a sparse matrix `A` using KLU[^ACM907]. -For sparse `A` with real or complex element type, the return type of `K` is -`KLUFactorization{Tv, Ti}`, with `Tv` = `Float64` or `ComplexF64` -respectively and `Ti` is an integer type (`Int32` or `Int64`). - +The values must be `Float64` or `ComplexF64`, and the indices `Int32` or `Int64`, as these are the +only combinations that the library provides. Anything else is a `MethodError`, so convert first, +for example with `klu(Float64.(A))`. The return type of `K` is `KLUFactorization{Tv, Ti}`. The individual components of the factorization `K` can be accessed by indexing: @@ -725,7 +722,8 @@ The relation between `K` and `A` is `K` further supports the following functions: -- `\\`, `ldiv!`, and [`KLU.solve!`](@ref) to solve linear systems +- `\\`, `ldiv!`, and [`KLU.solve!`](@ref) to solve linear systems, with right-hand sides of the same + element type as `K` - [`klu!`](@ref) to refactor with new values - [`rcond`](@ref), [`condest`](@ref), and [`rgrowth`](@ref) to check the quality of the factorization - `det`, `logabsdet`, and `logdet` @@ -740,15 +738,12 @@ The relation between `K` and `A` is Useful for cases where only the sparse structure of `A` is known at time of construction. !!! note - `klu(A::SparseMatrixCSC)` uses the KLU[^ACM907] library that is part of - SuiteSparse. As this library only supports sparse matrices with `Float64` or - `ComplexF64` elements, `klu` converts other values to one of these, in a new vector. + `klu(A::SparseMatrixCSC)` uses the KLU[^ACM907] library that is part of SuiteSparse. !!! warning "`K` shares vectors with the caller" To avoid using more memory, `K` does not copy what it is given. `klu(A)` keeps `nonzeros(A)`, - `klu(n, colptr, rowval, nzval)` keeps all three vectors, and `klu!(K, nzval)` keeps `nzval`, - unless the values had to be converted. It is the caller's responsibility not to change - them in a way that breaks `K`: + `klu(n, colptr, rowval, nzval)` keeps all three vectors, and `klu!(K, nzval)` keeps `nzval`. + It is the caller's responsibility not to change them in a way that breaks `K`: - Changing values does not change the current factors, but later calls to [`klu_factor!`](@ref), [`condest`](@ref), [`rgrowth`](@ref), and `copy` read the new values. @@ -773,25 +768,7 @@ function klu(n, check=true, allowsingular=false, full_factor=true, -) where {Ti<:KLUITypes, Tv<:AbstractFloat} - if Tv != Float64 - nzval = convert(Vector{Float64}, nzval) - end - K = KLUFactorization(n, colptr, rowval, nzval) - return full_factor ? klu_factor!(K; check, allowsingular) : klu_analyze!(K; check) -end - -function klu(n, - colptr::Vector{Ti}, - rowval::Vector{Ti}, - nzval::Vector{Tv}; - check=true, - allowsingular=false, - full_factor=true, -) where {Ti<:KLUITypes, Tv<:Complex} - if Tv != ComplexF64 - nzval = convert(Vector{ComplexF64}, nzval) - end +) where {Ti<:KLUITypes, Tv<:KLUTypes} K = KLUFactorization(n, colptr, rowval, nzval) return full_factor ? klu_factor!(K; check, allowsingular) : klu_analyze!(K; check) end @@ -800,7 +777,7 @@ function klu(A::SparseMatrixCSC{Tv, Ti}; check=true, allowsingular=false, full_factor=true, -) where {Tv<:Union{AbstractFloat, Complex}, Ti<:KLUITypes} +) where {Tv<:KLUTypes, Ti<:KLUITypes} n = size(A, 1) n == size(A, 2) || throw(DimensionMismatch()) return klu(n, @@ -825,9 +802,9 @@ If this throws because the new values are singular, `K` has no numeric factoriza later `klu!` or [`klu_factor!`](@ref) succeeds. Solving with `K` in the meantime factors the singular values again, and so throws as well. -For sparse `A` with real or complex element type, the return type of `K` is -`KLUFactorization{Tv, Ti}`, with `Tv` = `Float64` or `ComplexF64` -respectively and `Ti` is an integer type (`Int32` or `Int64`). +`A` must have the same value and index types as `K`, and `nzval` the same value type. Anything else +is a `MethodError`, so convert first. `K` keeps the vector it is given, see [`klu`](@ref) for what +this means. See also: [`klu`](@ref) @@ -837,14 +814,6 @@ See also: [`klu`](@ref) - `check::Bool`: If `true` (default) check for errors after the factorization. If `false` errors must be checked by the user with `klu.common.status`. - `allowsingular::Bool`: If `true` (default `false`) allow the factorization to proceed even if the matrix is singular. Note that this will allow for silent divide by zero errors in subsequent `solve!` or `ldiv!` calls if singularity is not checked by the user with `klu.common.status == KLU.KLU_SINGULAR` - -!!! note - `klu(A::SparseMatrixCSC)` uses the KLU[^ACM907] library that is part of - SuiteSparse. As this library only supports sparse matrices with `Float64` or - `ComplexF64` elements, `klu!` converts other values to one of these, in a new vector. - Otherwise `K` keeps the vector it is given, see [`klu`](@ref) for what this means. - -[^ACM907]: Davis, Timothy A., & Palamadai Natarajan, E. (2010). Algorithm 907: KLU, A Direct Sparse Solver for Circuit Simulation Problems. ACM Trans. Math. Softw., 37(3). doi:10.1145/1824801.1824814 """ klu! @@ -872,14 +841,6 @@ for Tv ∈ KLUValueTypes, Ti ∈ KLUIndexTypes end end -function klu!(K::AbstractKLUFactorization{ComplexF64}, nzval::Vector{U}; check=true, allowsingular=false) where {U<:Complex} - return klu!(K, convert(Vector{ComplexF64}, nzval); check, allowsingular) -end - -function klu!(K::AbstractKLUFactorization{Float64}, nzval::Vector{U}; check=true, allowsingular=false) where {U<:AbstractFloat} - return klu!(K, convert(Vector{Float64}, nzval); check, allowsingular) -end - # Compares without shifting the zero-based vector in place, as it may be shared with the # caller and with other factorizations, which could be inside libklu at the time. function _same_indices(zero_based::Vector, one_based::Vector) @@ -887,7 +848,7 @@ function _same_indices(zero_based::Vector, one_based::Vector) return all(i -> zero_based[i] + 1 == one_based[i], eachindex(zero_based, one_based)) end -function klu!(K::KLUFactorization{U}, S::SparseMatrixCSC{U}; check=true, allowsingular=false) where {U} +function klu!(K::KLUFactorization{Tv, Ti}, S::SparseMatrixCSC{Tv, Ti}; check=true, allowsingular=false) where {Tv, Ti} size(K) == size(S) || throw(ArgumentError("Sizes of K and S must match.")) # what should happen here when check = false? This is not really a KLU error code. _same_indices(K.colptr, S.colptr) && _same_indices(K.rowval, S.rowval) || @@ -966,7 +927,8 @@ end """ KLU.solve!(klu::Union{KLUFactorization, AdjointFactorization{Tv, KLUFactorization}, TransposeFactorization{Tv, KLUFactorization}}, B::StridedVecOrMat) -> B -Solve the linear system `AX = B` or `A'X = B` using the matrix factorization `klu` and right-hand side `B`. +Solve the linear system `AX = B` or `A'X = B` using the matrix factorization `klu` and right-hand side `B`, +which must have the same element type as `klu`. This function overwrites `B` with the solution `X`, for a new solution vector `X` use `klu \\ B`. @@ -1045,26 +1007,35 @@ for Tv ∈ KLUValueTypes, Ti ∈ KLUIndexTypes end end -function solve(klu, B; check=true) - X = copy(B) - return solve!(klu, X; check) -end - # we are not adding check to `ldiv!` since it is not in the contract. LinearAlgebra.ldiv!(klu::AbstractKLUFactorization{Tv}, B::StridedVecOrMat{Tv}) where {Tv<:KLUTypes} = solve!(klu, B) LinearAlgebra.ldiv!(klu::Union{AdjointFactorization{Tv, K},TransposeFactorization{Tv, K}}, B::StridedVecOrMat{Tv}) where {Tv, Ti, K<:AbstractKLUFactorization{Tv, Ti}} = solve!(klu, B) -function LinearAlgebra.ldiv!(klu::AbstractKLUFactorization{<:AbstractFloat}, B::StridedVecOrMat{<:Complex}) - imagX = solve(klu, imag(B)) - realX = solve(klu, real(B)) - map!(complex, B, realX, imagX) -end -function LinearAlgebra.ldiv!(klu::Union{AdjointFactorization{Tv, K},TransposeFactorization{Tv, K}}, B::StridedVecOrMat{<:Complex}) where {Tv<:AbstractFloat, Ti, K<:AbstractKLUFactorization{Tv, Ti}} - imagX = solve(klu, imag(B)) - realX = solve(klu, real(B)) - map!(complex, B, realX, imagX) +# A factorization, or one of the lazy wrappers that solve with it. +const _KLUSolver{Tv} = Union{ + AbstractKLUFactorization{Tv}, + AdjointFactorization{Tv, <:AbstractKLUFactorization{Tv}}, + TransposeFactorization{Tv, <:AbstractKLUFactorization{Tv}}, +} + +_rhs_error(::_KLUSolver{Tv}, ::Type{T}) where {Tv, T} = throw(ArgumentError( + "This factorization solves with right-hand sides of type $Tv, but got $T. Convert it first.")) + +# The generic methods in LinearAlgebra convert the right-hand side, which libklu does not provide. +(\)(K::_KLUSolver{Tv}, B::AbstractVecOrMat{Tv}) where {Tv} = ldiv!(K, copyto!(similar(B, Tv), B)) +(\)(K::_KLUSolver, B::AbstractVecOrMat{T}) where {T} = _rhs_error(K, T) +# LinearAlgebra has more specific methods for real factorizations with complex arrays. +for F in (:(AbstractKLUFactorization{T}), + :(AdjointFactorization{T, <:AbstractKLUFactorization{T}}), + :(TransposeFactorization{T, <:AbstractKLUFactorization{T}})) + @eval (\)(K::$F, B::VecOrMat{Complex{T}}) where {T<:LinearAlgebra.BlasReal} = _rhs_error(K, Complex{T}) end +LinearAlgebra.ldiv!(Y::AbstractVector{Tv}, K::_KLUSolver{Tv}, B::AbstractVector{Tv}) where {Tv} = ldiv!(K, copyto!(Y, B)) +LinearAlgebra.ldiv!(Y::AbstractMatrix{Tv}, K::_KLUSolver{Tv}, B::AbstractMatrix{Tv}) where {Tv} = ldiv!(K, copyto!(Y, B)) +LinearAlgebra.ldiv!(Y::AbstractVector{S}, K::_KLUSolver, B::AbstractVector{T}) where {S, T} = _rhs_error(K, S === T ? T : Union{S, T}) +LinearAlgebra.ldiv!(Y::AbstractMatrix{S}, K::_KLUSolver, B::AbstractMatrix{T}) where {S, T} = _rhs_error(K, S === T ? T : Union{S, T}) + end diff --git a/test/runtests.jl b/test/runtests.jl index e595dc9..99d69fc 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -29,12 +29,12 @@ using LinearAlgebra R = Diagonal(Tv == ComplexF64 ? complex.(klua.Rs) : klua.Rs) @test R \ A[klua.p, klua.q] ≈ (klua.L * klua.U + klua.F) - b = [8., 45., -3., 3., 19.] + b = Tv[8., 45., -3., 3., 19.] x = klua \ b @test x ≈ float([1:5;]) @test A*x ≈ b - z = complex.(b) + z = copy(b) x = ldiv!(klua, z) @test x ≈ float([1:5;]) @test z === x @@ -42,12 +42,12 @@ using LinearAlgebra # since klu_solve(, ) modifies , and has no field for . @test A*x ≈ b - b = [8., 20., 13., 6., 17.] + b = Tv[8., 20., 13., 6., 17.] x = klua'\b @test x ≈ float([1:5;]) @test A'*x ≈ b - z = complex.(b) + z = copy(b) x = ldiv!(adjoint(klua), z) @test x ≈ float([1:5;]) @test x === z @@ -57,17 +57,16 @@ using LinearAlgebra @test x ≈ float([1:5;]) @test transpose(A) * x ≈ b - x = ldiv!(transpose(klua), complex.(b)) + x = ldiv!(transpose(klua), copy(b)) @test x ≈ float([1:5;]) @test transpose(A) * x ≈ b - @inferred klua\fill(1, size(A, 2)) + @inferred klua\ones(Tv, size(A, 2)) for (op, M) in ((identity, A), (adjoint, A'), (transpose, transpose(A))) - X = complex.(rand(5, 2)) + X = rand(Tv, 5, 2) @test M * ldiv!(op(klua), copy(X)) ≈ X @test_throws DimensionMismatch ldiv!(op(klua), zeros(Tv, 6)) - @test_throws DimensionMismatch ldiv!(op(klua), zeros(ComplexF64, 6)) end @testset "Permutation vectors" begin @@ -355,13 +354,11 @@ end end @testset "full_factor = false" begin - for T in (Float64, ComplexF64, Float32, ComplexF32) + for T in (Float64, ComplexF64) A = sparse(Matrix{T}([1 0; 1 1])) nonzeros(A) .= zero(T) F = klu(A; full_factor = false) - # we do need both here--for non 64 bit types, the nzval of F is not the same as A's. - nonzeros(A) .= one(T) - nonzeros(F) .= one(T) + nonzeros(A) .= one(T) # F shares these values klu_factor!(F) b = ones(T, 2) @test F\b ≈ A\b @@ -570,7 +567,6 @@ end @test klu_factor!(K) \ b ≈ 2 * Ad \ b nonzeros(A) ./= 2 @test nonzeros(KLUFactorization(A)) !== nonzeros(A) - @test nonzeros(klu(SparseMatrixCSC{Float32, Int}(A))) isa Vector{Float64} colptr, rowval, nzval = [0, 1, 2, 4], [0, 1, 0, 2], [1.0, 2.0, 0.0, 3.0] K = klu(3, colptr, rowval, nzval) @@ -794,3 +790,71 @@ end @test D._symbolic == C_NULL && D._numeric == C_NULL @test K \ b ≈ x end + +@testset "Only the types that libklu provides" begin + A = sparse([2.0 1.0; 0.0 2.0]) + colptr, rowval = A.colptr .- 1, A.rowval .- 1 + K, Z = klu(A), klu(complex(A)) + + # Values + for T in (Float16, Float32, BigFloat, Int, Int32, Rational{Int}, ComplexF32, Complex{Int}, Complex{BigFloat}) + S = convert(SparseMatrixCSC{T, Int}, A) + v = nonzeros(S) + @test_throws MethodError klu(S) + @test_throws MethodError klu(2, colptr, rowval, v) + @test_throws MethodError KLUFactorization(S) + @test_throws MethodError KLUFactorization(2, colptr, rowval, v) + @test_throws MethodError klu!(K, S) + @test_throws MethodError klu!(Z, S) + @test_throws MethodError klu!(K, v) + @test_throws MethodError klu!(Z, v) + end + @test_throws MethodError klu(sparse([true false; false true])) + # A real and a complex factorization don't take each other's values. + @test_throws MethodError klu!(K, complex(A)) + @test_throws MethodError klu!(K, complex(nonzeros(A))) + @test_throws MethodError klu!(Z, A) + @test_throws MethodError klu!(Z, nonzeros(A)) + + # Indices, including a mix of the two that are provided. + @test_throws MethodError klu(convert(SparseMatrixCSC{Float64, Int16}, A)) + @test_throws MethodError klu(2, Int16.(colptr), Int16.(rowval), nonzeros(A)) + @test_throws MethodError klu(2, Int32.(colptr), rowval, nonzeros(A)) + @test_throws MethodError KLUFactorization(2, Int32.(colptr), rowval, nonzeros(A)) + @test_throws MethodError klu!(K, convert(SparseMatrixCSC{Float64, Int32}, A)) + + # All four combinations are accepted. + b = [1.0, 2.0] + for Tv in (Float64, ComplexF64), Ti in Base.uniontypes(KLUITypes) + S = convert(SparseMatrixCSC{Tv, Ti}, A) + F = klu(S) + @test F isa KLUFactorization{Tv, Ti} + @test klu!(F, S) \ Tv.(b) ≈ Matrix(A) \ b + @test inv(F) ≈ inv(Matrix(A)) + end + + # Right-hand sides + @testset "Right-hand sides with $op" for (op, M) in ((identity, A), (adjoint, A'), (transpose, transpose(A))) + for (F, Tv, others) in ((K, Float64, (Float16, Float32, Int, ComplexF32, ComplexF64)), + (Z, ComplexF64, (Float32, Float64, Int, ComplexF32, Complex{Int}))) + for T in others, B in (ones(T, 2), ones(T, 2, 3)) + @test_throws MethodError ldiv!(op(F), copy(B)) + @test_throws MethodError KLU.solve!(op(F), copy(B)) + @test_throws ArgumentError op(F) \ B + @test_throws ArgumentError ldiv!(copy(B), op(F), copy(B)) + @test_throws ArgumentError ldiv!(ones(Tv, size(B)), op(F), B) + @test_throws ArgumentError ldiv!(similar(B), op(F), ones(Tv, size(B))) + end + for B in (ones(Tv, 2), ones(Tv, 2, 3)) + @test M * (op(F) \ B) ≈ B + Y = similar(B) + @test ldiv!(Y, op(F), B) === Y + @test M * Y ≈ B + @test M * (op(F) \ view(B, :, axes(B)[2:end]...)) ≈ B + end + # Any array of the right element type can be given to `\`, which makes a copy. + V = view(Tv[1, 9, 2], [1, 3]) + @test M * (op(F) \ V) ≈ V + end + end +end