Merge pull request #32 from luke-kiernan/separate-sym-factor-rebase

Author: rayegun

docs/src/index.md

@@ -46,4 +46,8 @@ Please cite both [SuiteSparse](https://github.com/DrTimothyAldenDavis/SuiteSpars
 
 ```@docs
 klu
+klu!
+klu_refactor!
+nonzeros
+solve!
 ```

src/KLU.jl

@@ -2,9 +2,11 @@ module KLU
 
 using SparseArrays
 using SparseArrays: SparseMatrixCSC
-import SparseArrays: nnz
+import SparseArrays: nnz, nonzeros
 
 export klu, klu!
+export klu_factor!, klu_refactor!, solve!
+export nonzeros
 
 const libklu = :libklu
 include("wrappers.jl")
@@ -65,10 +67,10 @@ Data structure for parameters of and statistics generated by KLU functions.
 # Fields
 - `tol::Float64`: Partial pivoting tolerance for diagonal preference
 - `btf::Int64`: If `btf != 0` use BTF pre-ordering
-- `ordering::Int64`: If `ordering == 0` use AMD to permute, if `ordering == 1` use COLAMD,
+- `ordering::Int64`: If `ordering == 0` use AMD to permute, if `ordering == 1` use COLAMD, \
 if `ordering == 3` use the user provided ordering function.
-- `scale::Int64`: If `scale == 1` then `A[:,i] ./= sum(abs.(A[:,i]))`, if `scale == 2` then
-`A[:,i] ./= maximum(abs.(A[:,i]))`. If `scale == 0` no scaling is done, and the input is
+- `scale::Int64`: If `scale == 1` then `A[:,i] ./= sum(abs.(A[:,i]))`, if `scale == 2` then \
+`A[:,i] ./= maximum(abs.(A[:,i]))`. If `scale == 0` no scaling is done, and the input is \
 checked for errors if `scale >= 0`.
 
 See the [KLU User Guide](https://github.com/DrTimothyAldenDavis/SuiteSparse/raw/master/KLU/Doc/KLU_UserGuide.pdf)
@@ -185,6 +187,15 @@ function size(K::AbstractKLUFactorization, dim::Integer)
 end
 
 nnz(K::AbstractKLUFactorization) = K.lnz + K.unz + K.nzoff
+"""The nonzeros of the factorization `K`.
+!!! warning "nonzeros(K) may or may not point to nonzeros(A)"
+    After `K = klu(A)`, it may or may not be the case that `nonzeros(K) === nonzeros(A)`. \
+    If `A` is of type `SparseMatrixCSC{Float32, Int64}`, then `nonzeros(K)` will be \
+    a `Vector{Float64}`, while `nonzeros(A)` is a `Vector{Float32}`, so they're definitely distinct. \
+    On the other hand, if `A` has value type `Float64`, then changing the values in `nonzeros(A)` \
+    will change `nonzeros(K)` as well.
+"""
+nonzeros(K::AbstractKLUFactorization) = K.nzval
 
 if !isdefined(LinearAlgebra, :AdjointFactorization) # VERSION < v"1.10-"
     Base.adjoint(K::AbstractKLUFactorization) = Adjoint(K)
@@ -275,6 +286,7 @@ function getproperty(klu::AbstractKLUFactorization{Tv, Ti}, s::Symbol) where {Tv
         Lp = Vector{Ti}(undef, klu.n + 1)
         Li = Vector{Ti}(undef, lnz)
         Lx = Vector{Float64}(undef, lnz)
+        # could also be replaced with multiple dispatch.
         Lz = Tv == Float64 ? C_NULL : Vector{Float64}(undef, lnz)
         _extract!(klu; Lp, Li, Lx, Lz)
         return Lp, Li, Lx, Lz
@@ -375,6 +387,11 @@ function show(io::IO, mime::MIME{Symbol("text/plain")}, K::AbstractKLUFactorizat
     end
 end
 
+"""
+    klu_analyze!(K::KLUFactorization; check = true)
+
+Compute and store (as `K._symbolic`) a symbolic factorization of the sparse matrix \
+represented by the fields of `K`."""
 function klu_analyze!(K::KLUFactorization{Tv, Ti}; check=true) where {Tv, Ti<:KLUITypes}
     if K._symbolic != C_NULL return K end
     if Ti == Int64
@@ -390,7 +407,10 @@ function klu_analyze!(K::KLUFactorization{Tv, Ti}; check=true) where {Tv, Ti<:KL
     return K
 end
 
-# User provided permutation vectors:
+"""
+    klu_analyze!(K::KLUFactorization{Tv, Ti}, P::Vector{Ti}, Q::Vector{Ti}; check = true)
+
+Variant of `klu_analyze!` that allows for user-provided permutation permutation vectors `P` and `Q`."""
 function klu_analyze!(K::KLUFactorization{Tv, Ti}, P::Vector{Ti}, Q::Vector{Ti}; check=true) where {Tv, Ti<:KLUITypes}
     if K._symbolic != C_NULL return K end
     if Ti == Int64
@@ -496,11 +516,11 @@ existing `KLUFactorization` instance.
 # Arguments
   - `K::KLUFactorization`: The matrix factorization object to be factored.
   - `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`
+  - `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`
 
-The `K.common` struct can be used to modify certain options and parameters, see the
-[KLU documentation](https://github.com/DrTimothyAldenDavis/SuiteSparse/raw/master/KLU/Doc/KLU_UserGuide.pdf)
+The `K.common` struct can be used to modify certain options and parameters, see the \
+[KLU documentation](https://github.com/DrTimothyAldenDavis/SuiteSparse/raw/master/KLU/Doc/KLU_UserGuide.pdf) \
 or [`klu_common`](@ref) for more information.
 """
 klu_factor!
@@ -538,8 +558,10 @@ The relation between `K` and `A` is
 # Arguments
   - `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.
   - `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`
+  - `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`
+  - `full_factor::Bool`: if `true` (default), perform both numeric and symbolic factorization. If `false`, only perform symbolic factorization. \
+  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
@@ -549,31 +571,56 @@ The relation between `K` and `A` is
 
 [^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
 """
-function klu(n, colptr::Vector{Ti}, rowval::Vector{Ti}, nzval::Vector{Tv}; check=true, allowsingular=false) where {Ti<:KLUITypes, Tv<:AbstractFloat}
+function klu(n,
+    colptr::Vector{Ti},
+    rowval::Vector{Ti},
+    nzval::Vector{Tv};
+    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 klu_factor!(K; check, allowsingular)
+    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) where {Ti<:KLUITypes, Tv<:Complex}
+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
     K = KLUFactorization(n, colptr, rowval, nzval)
-    return klu_factor!(K; check, allowsingular)
+    return full_factor ? klu_factor!(K; check, allowsingular) : klu_analyze!(K; check)
 end
 
-function klu(A::SparseMatrixCSC{Tv, Ti}; check=true, allowsingular=false) where {Tv<:Union{AbstractFloat, Complex}, Ti<:KLUITypes}
+function klu(A::SparseMatrixCSC{Tv, Ti};
+    check=true,
+    allowsingular=false,
+    full_factor=true,
+) where {Tv<:Union{AbstractFloat, Complex}, Ti<:KLUITypes}
     n = size(A, 1)
     n == size(A, 2) || throw(DimensionMismatch())
-    return klu(n, decrement(A.colptr), decrement(A.rowval), A.nzval; check, allowsingular)
+    return klu(n,
+        decrement(A.colptr),
+        decrement(A.rowval),
+        A.nzval;
+        check,
+        allowsingular,
+        full_factor
+    )
 end
 
 """
     klu!(K::KLUFactorization, A::SparseMatrixCSC; check=true, allowsingular=false) -> K::KLUFactorization
-    klu(K::KLUFactorization, nzval::Vector{Tv}; check=true, allowsingular=false) -> K::KLUFactorization
+    klu!(K::KLUFactorization, nzval::Vector{Tv}; check=true, allowsingular=false) -> K::KLUFactorization
 
 Recompute the KLU factorization `K` using the values of `A` or `nzval`. The pattern of the original
 matrix used to create `K` must match the pattern of `A`. 
@@ -589,7 +636,7 @@ See also: [`klu`](@ref)
   - `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.
   - `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`
+    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
@@ -641,6 +688,12 @@ function klu!(K::KLUFactorization{U}, S::SparseMatrixCSC{U}; check=true, allowsi
     return klu!(K, S.nzval; check, allowsingular)
 end
 
+"""
+    klu_refactor!(...) -> K::KLUFactorization
+
+Same as [`klu!`](@ref): alias for naming consistency reasons."""
+klu_refactor!(args...) = klu!(args...)
+
 #B is the modified argument here. To match with the math it should be (klu, B). But convention is
 # modified first. Thoughts?
 
@@ -665,13 +718,13 @@ This function overwrites `B` with the solution `X`, for a new solution vector `X
 """
 solve!
 for Tv ∈ KLUValueTypes, Ti ∈ KLUIndexTypes
-    solve = _klu_name("solve", Tv, Ti)
+    solve_name = _klu_name("solve", Tv, Ti)
     @eval begin
         function solve!(klu::AbstractKLUFactorization{$Tv, $Ti}, B::StridedVecOrMat{$Tv}; check=true)
             stride(B, 1) == 1 || throw(ArgumentError("B must have unit strides"))
             klu._numeric == C_NULL && klu_factor!(klu)
             size(B, 1) == size(klu, 1) || throw(DimensionMismatch())
-            isok = $solve(klu._symbolic, klu._numeric, size(B, 1), size(B, 2), B, Ref(klu.common))
+            isok = $solve_name(klu._symbolic, klu._numeric, size(B, 1), size(B, 2), B, Ref(klu.common))
             isok == 0 && check && kluerror(klu.common)
             return B
         end

test/runtests.jl

@@ -2,7 +2,7 @@ using KLU
 using KLU: increment!, KLUITypes, decrement, klu!, KLUFactorization,
 klu_analyze!, klu_factor!
 using Test
-using SparseArrays: SparseMatrixCSC, sparse, nnz
+using SparseArrays: SparseMatrixCSC, sparse, nnz, nonzeros
 using LinearAlgebra
 
 @testset "KLU Wrappers" begin
@@ -123,3 +123,17 @@ end
         @test issuccess(klu(A; allowsingular=true); allowsingular=true)
     end
 end
+
+@testset "full_factor = false" begin
+    for T in (Float64, ComplexF64, Float32, ComplexF32)
+        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)
+        klu_factor!(F)
+        b = ones(T, 2)
+        @test F\b ≈ A\b
+    end
+end