Merge pull request #845 from hersle/check_pattern_first

Fix and speed up sparsity pattern check

Author: ChrisRackauckas

ext/LinearSolveSparseArraysExt.jl

@@ -201,7 +201,7 @@ function SciMLBase.solve!(
         cacheval = LinearSolve.@get_cacheval(cache, :UMFPACKFactorization)
         if alg.reuse_symbolic
             # Caches the symbolic factorization: https://github.com/JuliaLang/julia/pull/33738
-            if alg.check_pattern && pattern_changed(cacheval, A)
+            if length(cacheval.nzval) != length(A.nzval) || alg.check_pattern && pattern_changed(cacheval, A)
                 fact = lu(
                     SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),
                         nonzeros(A)),
@@ -331,7 +331,7 @@ function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::KLUFactorization;
     if cache.isfresh
         cacheval = LinearSolve.@get_cacheval(cache, :KLUFactorization)
         if alg.reuse_symbolic
-            if alg.check_pattern && pattern_changed(cacheval, A)
+            if length(cacheval.nzval) != length(A.nzval) || alg.check_pattern && pattern_changed(cacheval, A)
                 fact = KLU.klu(
                     SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),
                         nonzeros(A)),
@@ -455,9 +455,19 @@ function LinearSolve.pattern_changed(fact::Nothing, A::SparseArrays.SparseMatrix
 end
 
 function LinearSolve.pattern_changed(fact, A::SparseArrays.SparseMatrixCSC)
-    !(SparseArrays.decrement(SparseArrays.getcolptr(A)) ==
-      fact.colptr && SparseArrays.decrement(SparseArrays.getrowval(A)) ==
-      fact.rowval)
+    colptr0 = fact.colptr # has 0-based indices
+    colptr1 = SparseArrays.getcolptr(A) # has 1-based indices
+    length(colptr0) == length(colptr1) || return true
+    @inbounds for i in eachindex(colptr0)
+        colptr0[i] + 1 == colptr1[i] || return true
+    end
+    rowval0 = fact.rowval
+    rowval1 = SparseArrays.getrowval(A)
+    length(rowval0) == length(rowval1) || return true
+    @inbounds for i in eachindex(rowval0)
+        rowval0[i] + 1 == rowval1[i] || return true
+    end
+    return false
 end
 
 @static if Base.USE_GPL_LIBS

src/KLU/klu.jl

@@ -134,9 +134,9 @@ See the [`klu`](@ref) docs for more information.
 
 You typically should not construct this directly, instead use [`klu`](@ref).
 """
-mutable struct KLUFactorization{Tv <: KLUTypes, Ti <: KLUITypes} <:
+mutable struct KLUFactorization{Tv <: KLUTypes, Ti <: KLUITypes, Tklu <: Union{klu_l_common, klu_common}} <:
                AbstractKLUFactorization{Tv, Ti}
-    common::Union{klu_l_common, klu_common}
+    common::Tklu
     _symbolic::Ptr{Cvoid}
     _numeric::Ptr{Cvoid}
     n::Int
@@ -146,7 +146,7 @@ mutable struct KLUFactorization{Tv <: KLUTypes, Ti <: KLUITypes} <:
     function KLUFactorization(n, colptr, rowval, nzval)
         Ti = eltype(colptr)
         common = _common(Ti)
-        obj = new{eltype(nzval), Ti}(common, C_NULL, C_NULL, n, colptr, rowval, nzval)
+        obj = new{eltype(nzval), Ti, typeof(common)}(common, C_NULL, C_NULL, n, colptr, rowval, nzval)
         function f(klu)
             _free_symbolic(klu)
             _free_numeric(klu)