Backports for release 1.13 (#701)

- [x] #672 
- [x] #677 
- [x] #700 
- [x] #697 
- [x] #692 
- [x] #689 
- [x] #694 
- [x] #682

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

@@ -17,9 +17,9 @@ using LinearAlgebra: AdjOrTrans, AdjointFactorization, TransposeFactorization, m
 import Base: +, -, *, \, /, ==, zero
 import Base: Matrix, Vector
 import LinearAlgebra: mul!, ldiv!, rdiv!, cholesky, adjoint!, diag, eigen, dot,
-    issymmetric, istril, istriu, lu, tr, transpose!, tril!, triu!, isbanded,
+    issymmetric, istril, istriu, lu, tr, transpose!, tril!, triu!, isbanded, isdiag,
     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) ||
@@ -1742,7 +1745,7 @@ function opnormestinv(A::AbstractSparseMatrixCSC{T}, t::Integer = min(2,maximum(
                             repeated = true
                         end
                     end
-                    if !repeated
+                    if !repeated && 2^(n-1) ≥ 2t #we need enough non-parallel ±1 vectors
                         saux2 = S[1:n,j]' * S_old[1:n,1:t]
                         if _any_abs_eq(saux2,n)
                             repeated = true

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,86 +172,19 @@ function newcommon(; print = 0) # no printing from CHOLMOD by default
 end
 
 function getcommon(::Type{Int32})
-    init_once()
+    LibSuiteSparse.init_suitesparse()
     return get!(newcommon, task_local_storage(), :cholmod_common)::Ref{cholmod_common}
 end
 
 function getcommon(::Type{Int64})
-    init_once()
+    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)
 
-const init_once = Base.OncePerProcess{Nothing}() do
-    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 #
 ####################
@@ -1356,11 +1287,9 @@ function size(F::Factor, i::Integer)
     return 1
 end
 size(F::Factor) = (size(F, 1), size(F, 2))
-axes(A::Union{Dense,Sparse,Factor}) = map(Base.OneTo, size(A))
 
 IndexStyle(::Type{<:Dense}) = IndexLinear()
 
-size(FC::FactorComponent, i::Integer) = size(FC.F, i)
 size(FC::FactorComponent) = size(FC.F)
 
 adjoint(FC::FactorComponent{Tv,:L}) where {Tv} = FactorComponent{Tv,:U}(FC.F)

src/solvers/spqr.jl

@@ -115,7 +115,6 @@ struct QRSparseQ{Tv<:CHOLMOD.VTypes,Ti<:Integer} <: AbstractQ{Tv}
 end
 
 Base.size(Q::QRSparseQ) = (size(Q.factors, 1), size(Q.factors, 1))
-Base.axes(Q::QRSparseQ) = map(Base.OneTo, size(Q))
 
 Matrix{T}(Q::QRSparseQ) where {T} = lmul!(Q, Matrix{T}(I, size(Q, 1), min(size(Q, 1), Q.n)))
 
@@ -146,7 +145,6 @@ function Base.size(F::QRSparse, i::Integer)
         throw(ArgumentError("second argument must be positive"))
     end
 end
-Base.axes(F::QRSparse) = map(Base.OneTo, size(F))
 
 # From SPQR manual p. 6
 _default_tol(A::AbstractSparseMatrixCSC) =
@@ -510,8 +508,23 @@ function LinearAlgebra.ldiv!(X::StridedVecOrMat{T}, F::QRSparse{T}, B::StridedVe
         lmul!(adjoint(F.Q), W0)
 
         # Solve R*X = Q'*P*B
-        ldiv!(UpperTriangular(@view(F.R[Base.OneTo(rnk), Base.OneTo(rnk)])),
-              @view(W0[Base.OneTo(rnk), :]))
+        #
+        # We call generic_trimatdiv! directly instead of going through
+        # ldiv!(UpperTriangular(R_sub), ...) for two reasons:
+        # 1. UpperTriangular requires a square matrix, but F.R is m×n
+        #    so we can only take a column view R[:, 1:rnk] (which is
+        #    m×rnk, not square). A row+column view R[1:rnk, 1:rnk]
+        #    would be square but doesn't match SparseMatrixCSCView,
+        #    causing dispatch to a slow generic fallback.
+        # 2. generic_trimatdiv! is what UpperTriangular ldiv! dispatches
+        #    to anyway — calling it directly with uploc='U', isunitc='N',
+        #    tfun=identity is equivalent. The back-substitution loop
+        #    iterates over axes(B,1) = 1:rnk and searchsortedlast
+        #    excludes entries with row > j, so the extra rows in the
+        #    column view are never accessed.
+        W_rnk = @view(W0[Base.OneTo(rnk), :])
+        LinearAlgebra.generic_trimatdiv!(W_rnk, 'U', 'N', identity,
+                                         @view(F.R[:, Base.OneTo(rnk)]), W_rnk)
 
         # Apply right permutation: scatter solved rows into X using cpiv directly.
         # Zero X first so free variables (beyond rank) are zero in the basic solution.

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