Backports to release 1.12 (#702)

- [x] #689
- [x] #677 
- [x] #692

Author: dkarrasch

gen/Project.toml

@@ -4,5 +4,5 @@ JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
 SuiteSparse_jll = "bea87d4a-7f5b-5778-9afe-8cc45184846c"
 
 [compat]
-Clang = "0.18"
-JuliaFormatter = "1.0.45"
+Clang = "0.18, 0.19"
+JuliaFormatter = "2"

gen/generator.jl

@@ -22,6 +22,9 @@ else
 end
 include_dir = joinpath(artifact_dir, "include", "suitesparse") |> normpath
 
+config_h = joinpath(include_dir, "SuiteSparse_config.h")
+@assert isfile(config_h)
+
 cholmod_h = joinpath(include_dir, "cholmod.h")
 @assert isfile(cholmod_h)
 
@@ -42,7 +45,7 @@ options["general"]["output_file_path"] = joinpath(@__DIR__, "..", "src/solvers/w
 args = get_default_args()
 push!(args, "-I$include_dir")
 
-header_files = [cholmod_h, SuiteSparseQR_C_h, umfpack_h]
+header_files = [config_h, cholmod_h, SuiteSparseQR_C_h, umfpack_h]
 
 ctx = create_context(header_files, args, options)
 

src/SparseArrays.jl

@@ -19,7 +19,7 @@ import Base: Matrix, Vector
 import LinearAlgebra: mul!, ldiv!, rdiv!, cholesky, adjoint!, diag, eigen, dot,
     issymmetric, istril, istriu, lu, tr, transpose!, tril!, triu!, isbanded,
     cond, diagm, factorize, ishermitian, norm, opnorm, lmul!, rmul!, tril, triu,
-    matprod_dest, generic_matvecmul!, generic_matmatmul!, copytrito!
+    matprod_dest, generic_matvecmul!, generic_matmatmul!, generic_matmatmul_wrapper!, copytrito!
 
 import Base: adjoint, argmin, argmax, Array, broadcast, circshift!, complex, Complex,
     conj, conj!, convert, copy, copy!, copyto!, count, diff, findall, findmax, findmin,

src/linalg.jl

@@ -62,12 +62,12 @@ Base.@constprop :aggressive function spdensemul!(C, tA, tB, A, B, alpha, beta)
         _At_or_Ac_mul_B!(transpose, C, A, wrap(B, tB), alpha, beta)
     elseif tA_uc == 'C'
         _At_or_Ac_mul_B!(adjoint, C, A, wrap(B, tB), alpha, beta)
-    elseif tA_uc in ('S', 'H') && tB_uc == 'N'
+    elseif tA_uc in ('S', 'H')
         rangefun = isuppercase(tA) ? nzrangeup : nzrangelo
         diagop = tA_uc == 'S' ? identity : real
         odiagop = tA_uc == 'S' ? transpose : adjoint
         T = eltype(C)
-        _mul!(rangefun, diagop, odiagop, C, A, B, T(alpha), T(beta))
+        _mul!(rangefun, diagop, odiagop, C, A, wrap(B, tB), T(alpha), T(beta))
     else
         @stable_muladdmul LinearAlgebra._generic_matmatmul!(C, wrap(A, tA), wrap(B, tB), MulAddMul(alpha, beta))
     end
@@ -117,7 +117,7 @@ function _At_or_Ac_mul_B!(tfun::Function, C, A, B, α, β)
     C
 end
 
-Base.@constprop :aggressive function generic_matmatmul!(C::StridedMatrix, tA, tB, A::DenseMatrixUnion, B::SparseMatrixCSCUnion2, alpha::Number, beta::Number)
+Base.@constprop :aggressive function generic_matmatmul_wrapper!(C::StridedMatrix, tA, tB, A::DenseMatrixUnion, B::SparseMatrixCSCUnion2, alpha::Number, beta::Number, ::LinearAlgebra.BlasFlag.SyrkHerkGemm)
     transA = tA == 'N' ? identity : tA == 'T' ? transpose : adjoint
     if tB == 'N'
         _spmul!(C, transA(A), B, alpha, beta)
@@ -128,6 +128,9 @@ Base.@constprop :aggressive function generic_matmatmul!(C::StridedMatrix, tA, tB
     end
     return C
 end
+Base.@constprop :aggressive generic_matmatmul_wrapper!(C::StridedMatrix, tA, tB, A::DenseMatrixUnion, B::SparseMatrixCSCUnion2, alpha::Number, beta::Number, @nospecialize(val)) =
+    LinearAlgebra._generic_matmatmul!(C, wrap(A, tA), wrap(B, tB), alpha, beta)
+
 function _spmul!(C::StridedMatrix, X::DenseMatrixUnion, A::SparseMatrixCSCUnion2, α::Number, β::Number)
     mX, nX = size(X)
     nX == size(A, 1) ||

src/solvers/LibSuiteSparse.jl

@@ -1,12 +1,106 @@
 module LibSuiteSparse
 
 using SuiteSparse_jll
+import Libdl
 
 const TRUE  = Int32(1)
 const FALSE = Int32(0)
 
 include("wrappers.jl")
 
+const SUITESPARSE_MIN_VERSION = v"6.0.0"
+const BUILD_VERSION = VersionNumber(
+    SUITESPARSE_MAIN_VERSION,
+    SUITESPARSE_SUB_VERSION,
+    SUITESPARSE_SUBSUB_VERSION
+)
+
+public init_suitesparse
+
+"""
+    LibSuiteSparse.init_suitesparse
+
+Internal function which is used to initialize the SuiteSparse libraries to the correct
+    memory management functions. Any package which directly wraps one of the following
+    SuiteSparse libraries *must* ensure that this function is called before the use of that
+    library: AMD, CAMD, COLAMD, CCOLAMD, UMFPACK, CXSparse, CHOLMOD, KLU, BTF, LDL, RBio,
+    SPQR, SPEX, and ParU
+
+# Notes:
+- Currently this function only sets the memory management functions of SuiteSparse_config,
+    however there are also override functions for `printf`, `hypot`, and `divcomplex`.
+- SuiteSparse_config, and this initialization function, is not a dependency of CSparse,
+    GraphBLAS, or LAGraph.
+"""
+const init_suitesparse = Base.OncePerProcess{Nothing}() do
+    try
+        ### Check if the linked library is compatible with the Julia code
+        if Libdl.dlsym_e(Libdl.dlopen("libsuitesparseconfig"), :SuiteSparse_version) != C_NULL
+            current_version_array = Vector{Cint}(undef, 3)
+            SuiteSparse_version(current_version_array)
+            current_version = VersionNumber(current_version_array...)
+        else # SuiteSparse < 4.2.0 does not include SuiteSparse_version()
+            current_version = v"0.0.0"
+        end
+
+
+        if current_version < SUITESPARSE_MIN_VERSION
+            @warn """
+                SuiteSparse version incompatibility
+
+                Julia was compiled with SuiteSparse version $BUILD_VERSION. It is
+                currently linked with a version older than
+                $(SUITESPARSE_MIN_VERSION). This might cause Julia to
+                terminate when working with sparse matrix factorizations,
+                e.g. solving systems of equations with \\.
+
+                It is recommended that you use Julia with a recent version
+                of SuiteSparse, or download the generic binaries
+                from www.julialang.org, which ship with the correct
+                versions of all dependencies.
+                """
+        elseif BUILD_VERSION.major != current_version.major
+            @warn """
+                SuiteSparse version incompatibility
+
+                Julia was compiled with SuiteSparse version $BUILD_VERSION. It is
+                currently linked with version $current_version.
+                This might cause Julia to terminate when working with
+                sparse matrix factorizations, e.g. solving systems of
+                equations with \\.
+
+                It is recommended that you use Julia with the same major
+                version of SuiteSparse as the one used during the build, or
+                download the generic binaries from www.julialang.org,
+                which ship with the correct versions of all dependencies.
+                """
+        end
+
+        current_version >= v"6.0.0" && SuiteSparse_start()
+
+        # Register gc tracked allocator if SuiteSparse is new enough
+        if current_version >= v"7.0.0"
+            SuiteSparse_config_malloc_func_set(cglobal(:jl_malloc, Ptr{Cvoid}))
+            SuiteSparse_config_calloc_func_set(cglobal(:jl_calloc, Ptr{Cvoid}))
+            SuiteSparse_config_realloc_func_set(cglobal(:jl_realloc, Ptr{Cvoid}))
+            SuiteSparse_config_free_func_set(cglobal(:jl_free, Ptr{Cvoid}))
+        elseif current_version >= v"4.2.0"
+            cnfg = cglobal((:SuiteSparse_config, libsuitesparseconfig), Ptr{Cvoid})
+            unsafe_store!(cnfg, cglobal(:jl_malloc, Ptr{Cvoid}), 1)
+            unsafe_store!(cnfg, cglobal(:jl_calloc, Ptr{Cvoid}), 2)
+            unsafe_store!(cnfg, cglobal(:jl_realloc, Ptr{Cvoid}), 3)
+            unsafe_store!(cnfg, cglobal(:jl_free, Ptr{Cvoid}), 4)
+        end
+
+        current_version >= v"6.0.0" && atexit() do
+            SuiteSparse_finish()
+        end
+
+    catch ex
+        @error "Error during initialization of module LibSuiteSparse" exception=ex,catch_backtrace()
+    end
+end
+
 # exports
 const PREFIXES = ["cholmod_", "CHOLMOD_", "umfpack_"]
 for name in names(@__MODULE__; all=true), prefix in PREFIXES

src/solvers/cholmod.jl

@@ -24,8 +24,6 @@ import LinearAlgebra: (\), AdjointFactorization,
 
 using SparseArrays
 using SparseArrays: getcolptr, AbstractSparseVecOrMat
-import Libdl
-
 export
     Dense,
     Factor,
@@ -174,84 +172,19 @@ function newcommon(; print = 0) # no printing from CHOLMOD by default
 end
 
 function getcommon(::Type{Int32})
+    LibSuiteSparse.init_suitesparse()
     return get!(newcommon, task_local_storage(), :cholmod_common)::Ref{cholmod_common}
 end
 
 function getcommon(::Type{Int64})
+    LibSuiteSparse.init_suitesparse()
     return get!(newcommon_l, task_local_storage(), :cholmod_common_l)::Ref{cholmod_common}
 end
 
 getcommon() = getcommon(Int)
 
 const BUILD_VERSION = VersionNumber(CHOLMOD_MAIN_VERSION, CHOLMOD_SUB_VERSION, CHOLMOD_SUBSUB_VERSION)
 
-function __init__()
-    try
-        ### Check if the linked library is compatible with the Julia code
-        if Libdl.dlsym_e(Libdl.dlopen("libcholmod"), :cholmod_version) != C_NULL
-            current_version_array = Vector{Cint}(undef, 3)
-            cholmod_version(current_version_array)
-            current_version = VersionNumber(current_version_array...)
-        else # CHOLMOD < 2.1.1 does not include cholmod_version()
-            current_version = v"0.0.0"
-        end
-
-
-        if current_version < CHOLMOD_MIN_VERSION
-            @warn """
-                CHOLMOD version incompatibility
-
-                Julia was compiled with CHOLMOD version $BUILD_VERSION. It is
-                currently linked with a version older than
-                $(CHOLMOD_MIN_VERSION). This might cause Julia to
-                terminate when working with sparse matrix factorizations,
-                e.g. solving systems of equations with \\.
-
-                It is recommended that you use Julia with a recent version
-                of CHOLMOD, or download the generic binaries
-                from www.julialang.org, which ship with the correct
-                versions of all dependencies.
-                """
-        elseif BUILD_VERSION.major != current_version.major
-            @warn """
-                CHOLMOD version incompatibility
-
-                Julia was compiled with CHOLMOD version $BUILD_VERSION. It is
-                currently linked with version $current_version.
-                This might cause Julia to terminate when working with
-                sparse matrix factorizations, e.g. solving systems of
-                equations with \\.
-
-                It is recommended that you use Julia with the same major
-                version of CHOLMOD as the one used during the build, or
-                download the generic binaries from www.julialang.org,
-                which ship with the correct versions of all dependencies.
-                """
-        end
-
-        # Register gc tracked allocator if CHOLMOD is new enough
-        if current_version >= v"4.0.3"
-            ccall((:SuiteSparse_config_malloc_func_set, libsuitesparseconfig),
-                  Cvoid, (Ptr{Cvoid},), cglobal(:jl_malloc, Ptr{Cvoid}))
-            ccall((:SuiteSparse_config_calloc_func_set, libsuitesparseconfig),
-                  Cvoid, (Ptr{Cvoid},), cglobal(:jl_calloc, Ptr{Cvoid}))
-            ccall((:SuiteSparse_config_realloc_func_set, libsuitesparseconfig),
-                  Cvoid, (Ptr{Cvoid},), cglobal(:jl_realloc, Ptr{Cvoid}))
-            ccall((:SuiteSparse_config_free_func_set, libsuitesparseconfig),
-                  Cvoid, (Ptr{Cvoid},), cglobal(:jl_free, Ptr{Cvoid}))
-        elseif current_version >= v"3.0.0"
-            cnfg = cglobal((:SuiteSparse_config, libsuitesparseconfig), Ptr{Cvoid})
-            unsafe_store!(cnfg, cglobal(:jl_malloc, Ptr{Cvoid}), 1)
-            unsafe_store!(cnfg, cglobal(:jl_calloc, Ptr{Cvoid}), 2)
-            unsafe_store!(cnfg, cglobal(:jl_realloc, Ptr{Cvoid}), 3)
-            unsafe_store!(cnfg, cglobal(:jl_free, Ptr{Cvoid}), 4)
-        end
-
-    catch ex
-        @error "Error during initialization of module CHOLMOD" exception=ex,catch_backtrace()
-    end
-end
-
 ####################
 # Type definitions #
 ####################

src/solvers/umfpack.jl

@@ -325,15 +325,15 @@ is provided or `q` is `nothing`, UMFPACK's default is used. If the permutation i
 zero-based copy is made.
 
 The `control` vector defaults to the Julia SparseArrays package's default configuration for UMFPACK (NB: this is modified from the UMFPACK defaults to
-disable iterative refinement), but can be changed by passing a vector of length `UMFPACK_CONTROL`, see the UMFPACK manual for possible configurations. 
+disable iterative refinement), but can be changed by passing a vector of length `UMFPACK_CONTROL`, see the UMFPACK manual for possible configurations.
 For example to reenable iterative refinement:
 
     umfpack_control = SparseArrays.UMFPACK.get_umfpack_control(Float64, Int64) # read Julia default configuration for a Float64 sparse matrix
     SparseArrays.UMFPACK.show_umf_ctrl(umfpack_control) # optional - display values
     umfpack_control[SparseArrays.UMFPACK.JL_UMFPACK_IRSTEP] = 2.0 # reenable iterative refinement (2 is UMFPACK default max iterative refinement steps)
 
     Alu = lu(A; control = umfpack_control)
-    x = Alu \\ b   # solve Ax = b, including UMFPACK iterative refinement  
+    x = Alu \\ b   # solve Ax = b, including UMFPACK iterative refinement
 
 The individual components of the factorization `F` can be accessed by indexing:
 
@@ -1036,6 +1036,7 @@ for Tv in (:Float64, :ComplexF64), Ti in UmfpackIndexTypes
     # the control and info arrays
     _defaults = Symbol(umf_nm("defaults", Tv, Ti))
     @eval function get_umfpack_control(::Type{$Tv}, ::Type{$Ti})
+        LibSuiteSparse.init_suitesparse()
         control = Vector{Float64}(undef, UMFPACK_CONTROL)
         $_defaults(control)
         # Put julia's config here