Merge branch 'main' into export_abstractmatrixcsc

j-fu · j-fu · commit 8d5e1d4087f2 · 2024-12-05T11:29:02.000+01:00
diff --git a/.ci/test_and_change_uuid.jl b/.ci/test_and_change_uuid.jl
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
@@ -23,39 +23,27 @@ jobs:
         os:
           - ubuntu-latest
           - windows-latest
-          - macOS-13 # intel
         julia-arch:
           - x64
           - x86
-        exclude:
-          - os: macOS-13
-            julia-arch: x86
         include:
           - os: macOS-latest
             julia-arch: aarch64
             julia-version: 'nightly'
-            
+          - os: macOS-13
+            julia-arch: x64
+            julia-version: 'nightly'
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.julia-version }}
           arch: ${{ matrix.julia-arch }}
-      - uses: actions/cache@v4
-        env:
-          cache-name: cache-artifacts
-        with:
-          path: ~/.julia/artifacts
-          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
-          restore-keys: |
-            ${{ runner.os }}-test-${{ env.cache-name }}-
-            ${{ runner.os }}-test-${{ matrix.os }}
-            ${{ runner.os }}-
-      - run: julia --color=yes .ci/test_and_change_uuid.jl
+      - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v4
+      - uses: codecov/codecov-action@v5
         with:
           file: lcov.info
           token: ${{ secrets.CODECOV_TOKEN }}
@@ -80,7 +68,6 @@ jobs:
           version: ${{ matrix.julia-version }}
           arch: ${{ matrix.julia-arch }}
       - uses: julia-actions/cache@v2
-      - run: julia --color=yes .ci/test_and_change_uuid.jl
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
         env:
@@ -91,12 +78,10 @@ jobs:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2
         with:
-          # version: '1.6'
           version: 'nightly'
       - name: Generate docs
         run: |
-          julia --color=yes -e 'write("Project.toml", replace(read("Project.toml", String), r"uuid = .*?\n" =>"uuid = \"3f01184e-e22b-5df5-ae63-d93ebab69eaf\"\n"));'
-          julia --project --color=yes -e 'using Pkg; Pkg.activate("docs"); Pkg.instantiate(); Pkg.develop(PackageSpec(path = pwd()))'
+          julia --project --color=yes -e 'using Pkg; Pkg.activate("docs"); Pkg.develop(PackageSpec(path = pwd()))'
           julia --project=docs --color=yes docs/make.jl pdf
         env:
           DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}
diff --git a/LICENSE.md b/LICENSE.md
@@ -1,3 +1,22 @@
-The code in this repository is part of the Julia project. See the LICENSE file in the Julia project:
+MIT License
 
-https://github.com/JuliaLang/julia/blob/master/LICENSE.md
+Copyright (c) 2018-2024 SparseArrays.jl contributors:
+https://github.com/JuliaSparse/SparseArrays.jl/contributors
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/README.md b/README.md
@@ -20,16 +20,6 @@ This package ships as part of the Julia stdlib.
 
 SparseArrays.jl provides functionality for working with sparse arrays in Julia.
 
-## Using development versions of this package
-
-To use a newer version of this package, you need to build Julia from scratch. The build process is the same as any other build except that you need to change the commit used in `stdlib/SparseArrays.version`.
-
-It's also possible to load a development version of the package using [the trick used in the Section named "Using the development version of Pkg.jl" in the `Pkg.jl` repo](https://github.com/JuliaLang/Pkg.jl#using-the-development-version-of-pkgjl), but the capabilities are limited as all other packages will depend on the stdlib version of the package and will not work with the modified package. 
-
-The main environment may become inconsistent so you might need to run `Pkg.instantiate()` and/or `Pkg.resolve()` in the main or project environments if Julia complains about missing `Serialization.jl` in this package's dependencies. 
-
-For older (1.8 and before) `SuiteSparse.jl` needs to be bumped too.
-
 ## Updating SuiteSparse
 
 In order to upgrade SparseArrays.jl to use a new release of SuiteSparse, the following steps are necessary:
diff --git a/docs/src/index.md b/docs/src/index.md
@@ -206,20 +206,6 @@ section of the standard library reference.
 | [`sprandn(m,n,d)`](@ref)   | [`randn(m,n)`](@ref)   | Creates a *m*-by-*n* random matrix (of density *d*) with iid non-zero elements distributed according to the standard normal (Gaussian) distribution.                  |
 | [`sprandn(rng,m,n,d)`](@ref) | [`randn(rng,m,n)`](@ref) | Creates a *m*-by-*n* random matrix (of density *d*) with iid non-zero elements generated with the `rng` random number generator                                   |
 
-## [Sparse Linear Algebra](@id stdlib-sparse-linalg)
-
-Sparse matrix solvers call functions from [SuiteSparse](http://suitesparse.com). The following factorizations are available:
-
-1. [`cholesky`](@ref SparseArrays.CHOLMOD.cholesky)
-2. [`ldlt`](@ref SparseArrays.CHOLMOD.ldlt)
-3. [`lu`](@ref SparseArrays.UMFPACK.lu)
-4. [`qr`](@ref SparseArrays.SPQR.qr)
-
-| Type                  | Description                                   |
-|:----------------------|:--------------------------------------------- |
-| `CHOLMOD.Factor`      | Cholesky and LDLt factorizations              |
-| `UMFPACK.UmfpackLU`   | LU factorization                              |
-| `SPQR.QRSparse`       | QR factorization                              |
 
 ```@meta
 DocTestSetup = nothing
diff --git a/docs/src/solvers.md b/docs/src/solvers.md
@@ -4,6 +4,24 @@
 DocTestSetup = :(using LinearAlgebra, SparseArrays)
 ```
 
+## [Sparse Linear Algebra](@id stdlib-sparse-linalg)
+
+Sparse matrix solvers call functions from [SuiteSparse](http://suitesparse.com).
+
+The following factorizations are available:
+
+1. [`cholesky`](@ref SparseArrays.CHOLMOD.cholesky)
+2. [`ldlt`](@ref SparseArrays.CHOLMOD.ldlt)
+3. [`lu`](@ref SparseArrays.UMFPACK.lu)
+4. [`qr`](@ref SparseArrays.SPQR.qr)
+
+| Type                  | Description                                   |
+|:----------------------|:--------------------------------------------- |
+| `CHOLMOD.Factor`      | Cholesky and LDLt factorizations              |
+| `UMFPACK.UmfpackLU`   | LU factorization                              |
+| `SPQR.QRSparse`       | QR factorization                              |
+
+
 ```@docs; canonical=false
 SparseArrays.CHOLMOD.cholesky
 SparseArrays.CHOLMOD.cholesky!
diff --git a/src/abstractsparse.jl b/src/abstractsparse.jl
@@ -186,6 +186,13 @@ else
     end
 end
 
+macro allowscalar(p)
+    quote
+        $(allowscalar)($(esc(p)))
+        @Core.latestworld
+    end
+end
+
 @inline _is_fixed(::AbstractArray) = false
 @inline _is_fixed(A::AbstractArray, Bs::Vararg{Any,N}) where N = _is_fixed(A) || (N > 0 && _is_fixed(Bs...))
 macro if_move_fixed(a...)
diff --git a/src/linalg.jl b/src/linalg.jl
@@ -47,11 +47,11 @@ for op ∈ (:+, :-)
     end
 end
 
-generic_matmatmul!(C::StridedMatrix, tA, tB, A::SparseMatrixCSCUnion2, B::DenseMatrixUnion, alpha::Number, beta::Number) =
+@inline generic_matmatmul!(C::StridedMatrix, tA, tB, A::SparseMatrixCSCUnion2, B::DenseMatrixUnion, alpha::Number, beta::Number) =
     spdensemul!(C, tA, tB, A, B, alpha, beta)
-generic_matmatmul!(C::StridedMatrix, tA, tB, A::SparseMatrixCSCUnion2, B::AbstractTriangular, alpha::Number, beta::Number) =
+@inline generic_matmatmul!(C::StridedMatrix, tA, tB, A::SparseMatrixCSCUnion2, B::AbstractTriangular, alpha::Number, beta::Number) =
     spdensemul!(C, tA, tB, A, B, alpha, beta)
-generic_matvecmul!(C::StridedVecOrMat, tA, A::SparseMatrixCSCUnion2, B::DenseInputVector, alpha::Number, beta::Number) =
+@inline generic_matvecmul!(C::StridedVecOrMat, tA, A::SparseMatrixCSCUnion2, B::DenseInputVector, alpha::Number, beta::Number) =
     spdensemul!(C, tA, 'N', A, B, alpha, beta)
 
 Base.@constprop :aggressive function spdensemul!(C, tA, tB, A, B, alpha, beta)
@@ -188,6 +188,23 @@ function _A_mul_Bt_or_Bc!(tfun::Function, C::StridedMatrix, A::AbstractMatrix, B
     C
 end
 
+function *(A::Diagonal, b::AbstractSparseVector)
+    if size(A, 2) != length(b)
+        throw(
+            DimensionMismatch(lazy"The dimension of the matrix A $(size(A)) and of the vector b $(length(b))")
+        )
+    end
+    T = promote_eltype(A, b)
+    res = similar(b, T)
+    nzind_b = nonzeroinds(b)
+    nzval_b = nonzeros(b)
+    nzval_res = nonzeros(res)
+    for idx in eachindex(nzind_b)
+        nzval_res[idx] = A.diag[nzind_b[idx]] * nzval_b[idx]
+    end
+    return res
+end
+
 # Sparse matrix multiplication as described in [Gustavson, 1978]:
 # http://dl.acm.org/citation.cfm?id=355796
 
diff --git a/src/solvers/cholmod.jl b/src/solvers/cholmod.jl
@@ -901,6 +901,15 @@ Dense(A::StridedVecOrMatInclAdjAndTrans{T}) where
 
 Dense(A::Sparse) = sparse_to_dense(A)
 
+function Dense(ptr::Ptr{cholmod_dense})
+    if ptr == C_NULL
+        throw(ArgumentError("dense matrix construction failed for " *
+            "unknown reasons. Please submit a bug report."))
+    end
+    s = unsafe_load(ptr)
+    return Dense{jlxtype(s.xtype, s.dtype)}(ptr)
+end
+
 function Base.convert(::Type{Dense{Tnew}}, A::Dense{T}) where {Tnew, T}
     GC.@preserve A begin
         Ap = unsafe_load(pointer(A))
diff --git a/src/solvers/spqr.jl b/src/solvers/spqr.jl
@@ -204,10 +204,10 @@ function LinearAlgebra.qr(A::SparseMatrixCSC{Tv, Ti}; tol=_default_tol(A), order
         C_NULL, C_NULL, C_NULL, C_NULL,
         R, E, H, HPinv, HTau)
 
-    R_ = SparseMatrixCSC(Sparse(R[]))
-    return QRSparse(SparseMatrixCSC(Sparse(H[])),
-                    vec(Array(CHOLMOD.Dense{Tv}(HTau[]))),
-                    SparseMatrixCSC(min(size(A)...),
+    R_ = SparseMatrixCSC{Tv, Ti}(Sparse(R[]))
+    return QRSparse(SparseMatrixCSC{Tv, Ti}(Sparse(H[])),
+                    vec(Array{Tv}(CHOLMOD.Dense(HTau[]))),
+                    SparseMatrixCSC{Tv, Ti}(min(size(A)...),
                                     size(R_, 2),
                                     getcolptr(R_),
                                     rowvals(R_),
diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl
@@ -4167,7 +4167,7 @@ function is_hermsym(A::AbstractSparseMatrixCSC, check::Function)
     return true
 end
 
-function istriu(A::AbstractSparseMatrixCSC)
+function istriu(A::AbstractSparseMatrixCSC, k::Integer=0)
     m, n = size(A)
     colptr = getcolptr(A)
     rowval = rowvals(A)
@@ -4176,7 +4176,8 @@ function istriu(A::AbstractSparseMatrixCSC)
     for col = 1:min(n, m-1)
         l1 = colptr[col+1]-1
         for i = 0 : (l1 - colptr[col])
-            if rowval[l1-i] <= col
+            if rowval[l1-i] <= col - k
+                # rows preceeding the index would also lie above the band
                 break
             end
             if _isnotzero(nzval[l1-i])
@@ -4187,15 +4188,16 @@ function istriu(A::AbstractSparseMatrixCSC)
     return true
 end
 
-function istril(A::AbstractSparseMatrixCSC)
+function istril(A::AbstractSparseMatrixCSC, k::Integer=0)
     m, n = size(A)
     colptr = getcolptr(A)
     rowval = rowvals(A)
     nzval  = nonzeros(A)
 
     for col = 2:n
         for i = colptr[col] : (colptr[col+1]-1)
-            if rowval[i] >= col
+            if rowval[i] >= col - k
+                # subsequent rows would also lie below the band
                 break
             end
             if _isnotzero(nzval[i])
diff --git a/src/sparsevector.jl b/src/sparsevector.jl
@@ -1277,13 +1277,16 @@ function vcat(X1::_SparseConcatGroup, X::_SparseConcatGroup...)
     end
     return Base.typed_vcat(Base.promote_eltype(X1, X...), X1, X...)
 end
-function hvcat(rows::Tuple{Vararg{Int}}, X1::_SparseConcatGroup, X::_SparseConcatGroup...)
+function hvcat_internal(rows::Tuple{Vararg{Int}}, X1::_SparseConcatGroup, X::_SparseConcatGroup...)
     if anysparse(X1) || anysparse(X...)
         vcat(_hvcat_rows(rows, X1, X...)...)
     else
         Base.typed_hvcat(Base.promote_eltypeof(X1, X...), rows, X1, X...)
     end
 end
+function hvcat(rows::Tuple{Vararg{Int}}, X1::_SparseConcatGroup, X::_SparseConcatGroup...)
+    return hvcat_internal(rows, X1, X...)
+end
 function _hvcat_rows((row1, rows...)::Tuple{Vararg{Int}}, X::_SparseConcatGroup...)
     if row1 ≤ 0
         throw(ArgumentError("length of block row must be positive, got $row1"))
@@ -1304,9 +1307,8 @@ hcat(n1::Number, ns::Vararg{Number}) = invoke(hcat, Tuple{Vararg{Number}}, n1, n
 vcat(n1::Number, ns::Vararg{Number}) = invoke(vcat, Tuple{Vararg{Number}}, n1, ns...)
 hcat(n1::N, ns::Vararg{N}) where {N<:Number} = invoke(hcat, Tuple{Vararg{N}}, n1, ns...)
 vcat(n1::N, ns::Vararg{N}) where {N<:Number} = invoke(vcat, Tuple{Vararg{N}}, n1, ns...)
-hvcat(rows::Tuple{Vararg{Int}}, n1::Number, ns::Vararg{Number}) = invoke(hvcat, Tuple{typeof(rows), Vararg{Number}}, rows, n1, ns...)
-hvcat(rows::Tuple{Vararg{Int}}, n1::N, ns::Vararg{N}) where {N<:Number} = invoke(hvcat, Tuple{typeof(rows), Vararg{N}}, rows, n1, ns...)
-
+hvcat(rows::Tuple{Vararg{Int}}, n1::Number, ns::Vararg{Number}) = hvcat_internal(rows, n1, ns...)
+hvcat(rows::Tuple{Vararg{Int}}, n1::N, ns::Vararg{N}) where {N<:Number} = hvcat_internal(rows, n1, ns...)
 
 ### Efficient repetition of sparse vectors
 
diff --git a/test/allowscalar.jl b/test/allowscalar.jl
@@ -4,20 +4,20 @@ using Test, SparseArrays
     A = sprandn(10, 20, 0.9)
     A[1, 1] = 2
     @test A[1, 1] == 2
-    SparseArrays.allowscalar(false)
+    SparseArrays.@allowscalar(false)
     @test_throws Any A[1, 1]
     @test_throws Any A[1, 1] = 2
-    SparseArrays.allowscalar(true)
+    SparseArrays.@allowscalar(true)
     @test A[1, 1] == 2
     A[1, 1] = 3
     @test A[1, 1] == 3
 
     B = sprandn(10, 0.9)
     B[1] = 2
     @test B[1] == 2
-    SparseArrays.allowscalar(false)
+    SparseArrays.@allowscalar(false)
     @test_throws Any B[1]
-    SparseArrays.allowscalar(true)
+    SparseArrays.@allowscalar(true)
     @test B[1] == 2
     B[1] = 3
     @test B[1] == 3
diff --git a/test/linalg.jl b/test/linalg.jl
@@ -673,6 +673,21 @@ end
     end
 end
 
+@testset "diagonal - sparse vector mutliplication" begin
+    for _ in 1:10
+        b = spzeros(10)
+        b[1:3] .= 1:3
+        A = Diagonal(randn(10))
+        @test norm(A * b - A * Vector(b)) <= 10eps()
+        @test norm(A * b - Array(A) * b) <= 10eps()
+        Ac = Diagonal(randn(Complex{Float64}, 10))
+        @test norm(Ac * b - Ac * Vector(b)) <= 10eps()
+        @test norm(Ac * b - Array(Ac) * b) <= 10eps()
+        @test_throws DimensionMismatch A * [b; 1]
+        @test_throws DimensionMismatch A * b[1:end-1]
+    end
+end
+
 @testset "sparse matrix * BitArray" begin
     A = sprand(5,5,0.3)
     MA = Array(A)
diff --git a/test/sparsematrix_ops.jl b/test/sparsematrix_ops.jl
diff --git a/test/spqr.jl b/test/spqr.jl
