Skip to content

Commit 0d1ecda

Browse files
rayegunrayegun
committed
Actually commit before pushing :)
1 parent d9fcc3d commit 0d1ecda

File tree

2 files changed

+122
-36
lines changed

2 files changed

+122
-36
lines changed

‎src/KLU.jl

Lines changed: 113 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -355,8 +355,8 @@ function getproperty(klu::AbstractKLUFactorization{Tv, Ti}, s::Symbol) where {Tv
355355
end
356356
end
357357

358-
function LinearAlgebra.issuccess(K::AbstractKLUFactorization)
359-
return K.common.status == KLU_OK && K._numeric != C_NULL
358+
function LinearAlgebra.issuccess(K::AbstractKLUFactorization; allowsingular=false)
359+
return (allowsingular ? K.common.status >= KLU_OK : K.common.status == KLU_OK) && K._numeric != C_NULL
360360
end
361361
function show(io::IO, mime::MIME{Symbol("text/plain")}, K::AbstractKLUFactorization)
362362
summary(io, K); println(io)
@@ -375,14 +375,14 @@ function show(io::IO, mime::MIME{Symbol("text/plain")}, K::AbstractKLUFactorizat
375375
end
376376
end
377377

378-
function klu_analyze!(K::KLUFactorization{Tv, Ti}) where {Tv, Ti<:KLUITypes}
378+
function klu_analyze!(K::KLUFactorization{Tv, Ti}; check=true) where {Tv, Ti<:KLUITypes}
379379
if K._symbolic != C_NULL return K end
380380
if Ti == Int64
381381
sym = klu_l_analyze(K.n, K.colptr, K.rowval, Ref(K.common))
382382
else
383383
sym = klu_analyze(K.n, K.colptr, K.rowval, Ref(K.common))
384384
end
385-
if sym == C_NULL
385+
if sym == C_NULL && check
386386
kluerror(K.common)
387387
else
388388
K._symbolic = sym
@@ -391,14 +391,14 @@ function klu_analyze!(K::KLUFactorization{Tv, Ti}) where {Tv, Ti<:KLUITypes}
391391
end
392392

393393
# User provided permutation vectors:
394-
function klu_analyze!(K::KLUFactorization{Tv, Ti}, P::Vector{Ti}, Q::Vector{Ti}) where {Tv, Ti<:KLUITypes}
394+
function klu_analyze!(K::KLUFactorization{Tv, Ti}, P::Vector{Ti}, Q::Vector{Ti}; check=true) where {Tv, Ti<:KLUITypes}
395395
if K._symbolic != C_NULL return K end
396396
if Ti == Int64
397397
sym = klu_l_analyze_given(K.n, K.colptr, K.rowval, P, Q, Ref(K.common))
398398
else
399399
sym = klu_analyze_given(K.n, K.colptr, K.rowval, P, Q, Ref(K.common))
400400
end
401-
if sym == C_NULL
401+
if sym == C_NULL && check
402402
kluerror(K.common)
403403
else
404404
K._symbolic = sym
@@ -409,14 +409,25 @@ end
409409
for Tv ∈ KLUValueTypes, Ti ∈ KLUIndexTypes
410410
factor = _klu_name("factor", Tv, Ti)
411411
@eval begin
412-
function klu_factor!(K::KLUFactorization{$Tv, $Ti})
413-
K._symbolic == C_NULL && klu_analyze!(K)
414-
num = $factor(K.colptr, K.rowval, K.nzval, K._symbolic, Ref(K.common))
415-
if num == C_NULL
412+
function klu_factor!(K::KLUFactorization{$Tv, $Ti}; check=true, allowsingular=false)
413+
K._symbolic == C_NULL && K.common.status >= KLU_OK && klu_analyze!(K)
414+
if K._symbolic != C_NULL && K.common.status >= KLU_OK
415+
K.common.halt_if_singular = !allowsingular && check
416+
num = $factor(K.colptr, K.rowval, K.nzval, K._symbolic, Ref(K.common))
417+
K.common.halt_if_singular = true
418+
else
419+
num = C_NULL
420+
end
421+
if num == C_NULL && check
416422
kluerror(K.common)
417423
else
418-
K._numeric = num
424+
if allowsingular
425+
K.common.status < KLU_OK && check && kluerror(K.common)
426+
else
427+
(K.common.status == KLU_OK) || (check && kluerror(K.common))
428+
end
419429
end
430+
K._numeric = num
420431
return K
421432
end
422433
end
@@ -476,12 +487,18 @@ end
476487

477488

478489
"""
479-
klu_factor!(K::KLUFactorization)
490+
klu_factor!(K::KLUFactorization; check=true, allowsingular=false) -> K::KLUFactorization
480491
481492
Factor `K` into components `K.L`, `K.U`, and `K.F`.
482493
This function will perform both the symbolic and numeric steps of factoriation on an
483494
existing `KLUFactorization` instance.
484495
496+
# Arguments
497+
- `K::KLUFactorization`: The matrix factorization object to be factored.
498+
- `check::Bool`: If `true` (default) check for errors after the factorization. If `false` errors must be checked by the user with `klu.common.status`.
499+
- `allowsingular::Bool`: If `true` (default `false`) allow the factorization to proceed even if the matrix is singular. Note that this will allow for
500+
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`
501+
485502
The `K.common` struct can be used to modify certain options and parameters, see the
486503
[KLU documentation](https://github.com/DrTimothyAldenDavis/SuiteSparse/raw/master/KLU/Doc/KLU_UserGuide.pdf)
487504
or [`klu_common`](@ref) for more information.
@@ -518,6 +535,12 @@ The relation between `K` and `A` is
518535
519536
- `LinearAlgebra.\\`
520537
538+
# Arguments
539+
- `A::SparseMatrixCSC` or `n::Integer`, `colptr::Vector{Ti}`, `rowval::Vector{Ti}`, `nzval::Vector{Tv}`: The sparse matrix or the zero-based sparse matrix components to be factored.
540+
- `check::Bool`: If `true` (default) check for errors after the factorization. If `false` errors must be checked by the user with `klu.common.status`.
541+
- `allowsingular::Bool`: If `true` (default `false`) allow the factorization to proceed even if the matrix is singular. Note that this will allow for
542+
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`
543+
521544
!!! note
522545
`klu(A::SparseMatrixCSC)` uses the KLU[^ACM907] library that is part of
523546
SuiteSparse. As this library only supports sparse matrices with [`Float64`](@ref) or
@@ -526,36 +549,68 @@ The relation between `K` and `A` is
526549
527550
[^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
528551
"""
529-
function klu(n, colptr::Vector{Ti}, rowval::Vector{Ti}, nzval::Vector{Tv}) where {Ti<:KLUITypes, Tv<:AbstractFloat}
552+
function klu(n, colptr::Vector{Ti}, rowval::Vector{Ti}, nzval::Vector{Tv}; check=true, allowsingular=false) where {Ti<:KLUITypes, Tv<:AbstractFloat}
530553
if Tv != Float64
531554
nzval = convert(Vector{Float64}, nzval)
532555
end
533556
K = KLUFactorization(n, colptr, rowval, nzval)
534-
return klu_factor!(K)
557+
return klu_factor!(K; check, allowsingular)
535558
end
536559

537-
function klu(n, colptr::Vector{Ti}, rowval::Vector{Ti}, nzval::Vector{Tv}) where {Ti<:KLUITypes, Tv<:Complex}
560+
function klu(n, colptr::Vector{Ti}, rowval::Vector{Ti}, nzval::Vector{Tv}; check=true, allowsingular=false) where {Ti<:KLUITypes, Tv<:Complex}
538561
if Tv != ComplexF64
539562
nzval = convert(Vector{ComplexF64}, nzval)
540563
end
541564
K = KLUFactorization(n, colptr, rowval, nzval)
542-
return klu_factor!(K)
565+
return klu_factor!(K; check, allowsingular)
543566
end
544567

545-
function klu(A::SparseMatrixCSC{Tv, Ti}) where {Tv<:Union{AbstractFloat, Complex}, Ti<:KLUITypes}
568+
function klu(A::SparseMatrixCSC{Tv, Ti}; check=true, allowsingular=false) where {Tv<:Union{AbstractFloat, Complex}, Ti<:KLUITypes}
546569
n = size(A, 1)
547570
n == size(A, 2) || throw(DimensionMismatch())
548-
return klu(n, decrement(A.colptr), decrement(A.rowval), A.nzval)
571+
return klu(n, decrement(A.colptr), decrement(A.rowval), A.nzval; check, allowsingular)
549572
end
550573

574+
"""
575+
klu!(K::KLUFactorization, A::SparseMatrixCSC; check=true, allowsingular=false) -> K::KLUFactorization
576+
klu(K::KLUFactorization, nzval::Vector{Tv}; check=true, allowsingular=false) -> K::KLUFactorization
577+
578+
Recompute the KLU factorization `K` using the values of `A` or `nzval`. The pattern of the original
579+
matrix used to create `K` must match the pattern of `A`.
580+
581+
For sparse `A` with real or complex element type, the return type of `K` is
582+
`KLUFactorization{Tv, Ti}`, with `Tv` = `Float64` or `ComplexF64`
583+
respectively and `Ti` is an integer type (`Int32` or `Int64`).
584+
585+
See also: [`klu`](@ref)
586+
587+
# Arguments
588+
- `K::KLUFactorization`: The matrix factorization object, previously created by a call to `klu`, to be re-factored.
589+
- `A::SparseMatrixCSC` or `n::Integer`, `colptr::Vector{Ti}`, `rowval::Vector{Ti}`, `nzval::Vector{Tv}`: The sparse matrix or the zero-based sparse matrix components to be factored.
590+
- `check::Bool`: If `true` (default) check for errors after the factorization. If `false` errors must be checked by the user with `klu.common.status`.
591+
- `allowsingular::Bool`: If `true` (default `false`) allow the factorization to proceed even if the matrix is singular. Note that this will allow for
592+
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`
593+
594+
!!! note
595+
`klu(A::SparseMatrixCSC)` uses the KLU[^ACM907] library that is part of
596+
SuiteSparse. As this library only supports sparse matrices with [`Float64`](@ref) or
597+
`ComplexF64` elements, `lu` converts `A` into a copy that is of type
598+
`SparseMatrixCSC{Float64}` or `SparseMatrixCSC{ComplexF64}` as appropriate.
599+
600+
[^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
601+
"""
602+
klu!
603+
551604
for Tv ∈ KLUValueTypes, Ti ∈ KLUIndexTypes
552605
refactor = _klu_name("refactor", Tv, Ti)
553606
@eval begin
554-
function klu!(K::KLUFactorization{$Tv, $Ti}, nzval::Vector{$Tv})
607+
function klu!(K::KLUFactorization{$Tv, $Ti}, nzval::Vector{$Tv}; check=true, allowsingular=false)
555608
length(nzval) != length(K.nzval) && throw(DimensionMismatch())
556609
K.nzval = nzval
610+
K.common.halt_if_singular = !allowsingular && check
557611
ok = $refactor(K.colptr, K.rowval, K.nzval, K._symbolic, K._numeric, Ref(K.common))
558-
if ok == 1
612+
K.common.halt_if_singular = true
613+
if (ok == 1 || !check || (allowsingular && K.common.status >= KLU_OK))
559614
return K
560615
else
561616
kluerror(K.common)
@@ -564,39 +619,60 @@ for Tv ∈ KLUValueTypes, Ti ∈ KLUIndexTypes
564619
end
565620
end
566621

567-
function klu!(K::AbstractKLUFactorization{ComplexF64}, nzval::Vector{U}) where {U<:Complex}
568-
return klu!(K, convert(Vector{ComplexF64}, nzval))
622+
function klu!(K::AbstractKLUFactorization{ComplexF64}, nzval::Vector{U}; check=true, allowsingular=false) where {U<:Complex}
623+
return klu!(K, convert(Vector{ComplexF64}, nzval); check, allowsingular)
569624
end
570625

571-
function klu!(K::AbstractKLUFactorization{Float64}, nzval::Vector{U}) where {U<:AbstractFloat}
572-
return klu!(K, convert(Vector{Float64}, nzval))
626+
function klu!(K::AbstractKLUFactorization{Float64}, nzval::Vector{U}; check=true, allowsingular=false) where {U<:AbstractFloat}
627+
return klu!(K, convert(Vector{Float64}, nzval); check, allowsingular)
573628
end
574629

575-
function klu!(K::KLUFactorization{U}, S::SparseMatrixCSC{U}) where {U}
630+
function klu!(K::KLUFactorization{U}, S::SparseMatrixCSC{U}; check=true, allowsingular=false) where {U}
576631
size(K) == size(S) || throw(ArgumentError("Sizes of K and S must match."))
577632
increment!(K.colptr)
578633
increment!(K.rowval)
634+
# what should happen here when check = false? This is not really a KLU error code.
579635
K.colptr == S.colptr && K.rowval == S.rowval ||
580636
(decrement!(K.colptr); decrement!(K.rowval);
581637
throw(ArgumentError("The pattern of the original matrix must match the pattern of the refactor."))
582638
)
583639
decrement!(K.colptr)
584640
decrement!(K.rowval)
585-
return klu!(K, S.nzval)
641+
return klu!(K, S.nzval; check, allowsingular)
586642
end
587643

588-
589644
#B is the modified argument here. To match with the math it should be (klu, B). But convention is
590645
# modified first. Thoughts?
646+
647+
# with check=false it *is technically* possible to enter a seriously invalid state wherein `klu.common`
648+
# is undefined. I don't think this is possible in practice though, since we throw on invalid `common` on construction.
649+
"""
650+
solve!(klu::Union{KLUFactorization, AdjointFact{Tv, KLUFactorization}, TransposeFact{Tv, KLUFactorization}}, B::StridedVecOrMat) -> B
651+
652+
Solve the linear system `AX = B` or `A'X = B` where using the matrix facotrization `klu` and right-hand side `B`.
653+
654+
This function overwrites `B` with the solution `X`, for a new solution vector `X` use `solve(klu, B)`.
655+
656+
# Arguments
657+
- `klu::KLUFactorization`: The matrix factorization of `A` to use in the solution.
658+
- `B::StridedVecOrMat`: The right-hand side of the linear system.
659+
- `check::Bool`: If `true` (default) check for errors after the solve. If `false` errors must be checked using `klu.common.status`. The return value of this function is always `B`, so the status of the solve *must* be checked using the factorization object itself when `check=false`.
660+
661+
!!! warning "Solve with a singular factorization object"
662+
If the factorization object `klu` has `klu.common.status == KLU.KLU_SINGULAR` then the `solve!` or `ldiv!` will result in a silent divide by zero error.
663+
664+
This status should be checked by the user before solve calls if singularity checks were disabled on factorization using `check=false` or `allowsingular=true`.
665+
"""
666+
solve!
591667
for Tv ∈ KLUValueTypes, Ti ∈ KLUIndexTypes
592668
solve = _klu_name("solve", Tv, Ti)
593669
@eval begin
594-
function solve!(klu::AbstractKLUFactorization{$Tv, $Ti}, B::StridedVecOrMat{$Tv})
670+
function solve!(klu::AbstractKLUFactorization{$Tv, $Ti}, B::StridedVecOrMat{$Tv}; check=true)
595671
stride(B, 1) == 1 || throw(ArgumentError("B must have unit strides"))
596672
klu._numeric == C_NULL && klu_factor!(klu)
597673
size(B, 1) == size(klu, 1) || throw(DimensionMismatch())
598674
isok = $solve(klu._symbolic, klu._numeric, size(B, 1), size(B, 2), B, Ref(klu.common))
599-
isok == 0 && kluerror(klu.common)
675+
isok == 0 && check && kluerror(klu.common)
600676
return B
601677
end
602678
end
@@ -610,33 +686,35 @@ for Tv ∈ KLUValueTypes, Ti ∈ KLUIndexTypes
610686
call = :($tsolve(klu._symbolic, klu._numeric, size(B, 1), size(B, 2), B, Ref(klu.common)))
611687
end
612688
@eval begin
613-
function solve!(klu::AdjointFact{$Tv, K}, B::StridedVecOrMat{$Tv}) where {K<:AbstractKLUFactorization{$Tv, $Ti}}
689+
function solve!(klu::AdjointFact{$Tv, K}, B::StridedVecOrMat{$Tv}; check=true) where {K<:AbstractKLUFactorization{$Tv, $Ti}}
614690
conj = 1
615691
klu = parent(klu)
616692
stride(B, 1) == 1 || throw(ArgumentError("B must have unit strides"))
617693
klu._numeric == C_NULL && klu_factor!(klu)
618694
size(B, 1) == size(klu, 1) || throw(DimensionMismatch())
619695
isok = $call
620-
isok == 0 && kluerror(klu.common)
696+
isok == 0 && check && kluerror(klu.common)
621697
return B
622698
end
623-
function solve!(klu::TransposeFact{$Tv, K}, B::StridedVecOrMat{$Tv}) where {K<: AbstractKLUFactorization{$Tv, $Ti}}
699+
function solve!(klu::TransposeFact{$Tv, K}, B::StridedVecOrMat{$Tv}; check=true) where {K<: AbstractKLUFactorization{$Tv, $Ti}}
624700
conj = 0
625701
klu = parent(klu)
626702
stride(B, 1) == 1 || throw(ArgumentError("B must have unit strides"))
627703
klu._numeric == C_NULL && klu_factor!(klu)
628704
size(B, 1) == size(klu, 1) || throw(DimensionMismatch())
629705
isok = $call
630-
isok == 0 && kluerror(klu.common)
706+
isok == 0 && check && kluerror(klu.common)
631707
return B
632708
end
633709
end
634710
end
635711

636-
function solve(klu, B)
712+
function solve(klu, B; check=true)
637713
X = copy(B)
638-
return solve!(klu, X)
714+
return solve!(klu, X; check)
639715
end
716+
717+
# we are not adding check to `ldiv!` since it is not in the contract.
640718
LinearAlgebra.ldiv!(klu::AbstractKLUFactorization{Tv}, B::StridedVecOrMat{Tv}) where {Tv<:KLUTypes} =
641719
solve!(klu, B)
642720
LinearAlgebra.ldiv!(klu::Union{AdjointFact{Tv, K},TransposeFact{Tv, K}}, B::StridedVecOrMat{Tv}) where {Tv, Ti, K<:AbstractKLUFactorization{Tv, Ti}} =

‎test/runtests.jl

Lines changed: 9 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -114,4 +114,12 @@ end
114114
[0.2775474841561938, 0.19549953706849155, 0.7221976371086005, 0.4339373082200655, 0.983079431343046, 0.10918778088879799, 0.3728676112188065, 0.9134045061777432, 0.14560891622463457, 0.7715210431553383, 0.2945501295372417, 0.6134722502122134, 0.8777181195348973, 0.6995382425541914, 0.9562490972786235, 0.27001502642215325, 0.8573661029146233, 0.13020432565448115, 0.9221068910751316, 0.17087414970038983, 0.7062975193151109, 0.7668596005709167, 0.46967704631299334, 0.31764226298560083, 0.39054386892157833, 0.36610203401046015, 0.16689896372140534, 0.9624322297755521, 0.1478381603984824, 0.45423514524961806, 0.5564610482579242, 0.1844671322175948, 0.0823893170743808, 0.25409993152712307, 0.10475245199943273, 0.5863595004922162, 0.14733912690513562, 0.6504152422320895, 0.4339054908933866, 0.27614384058497166, 0.4019619228414033, 0.8631491210976487, 0.20159747073826084, 0.3273367915690062, 0.23866880928640288, 0.9557759456784265, 0.016351125161178537, 0.5320355909884844, 0.9010930260468242, 0.3780686420068593, 0.6375477164214856, 0.9850645305956225, 0.5366242762582065, 0.08835652070698918, 0.9877090693717305, 0.9775298646022268, 0.9759511830494418]
115115
)
116116
@test klu(A).F.nzval ≈ [0.20623152093700403, 0.09038754099133628, 1.0]
117-
end
117+
end
118+
119+
@testset "check=false" begin
120+
for A in sparse.((Float64[1 2; 0 0], ComplexF64[1 2; 0 0]))
121+
@test_throws SingularException klu(A)
122+
@test !issuccess(klu(A; check = false))
123+
@test issuccess(klu(A; allowsingular=true); allowsingular=true)
124+
end
125+
end

0 commit comments

Comments
 (0)