diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 232c421d..9719cb6e 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -30,11 +30,11 @@ jobs:
           - os: macOS-latest
             julia-arch: aarch64
             julia-version: 'nightly'
-          - os: macOS-13
+          - os: macOS-15-intel
             julia-arch: x64
             julia-version: 'nightly'
     steps:
-      - uses: actions/checkout@v5
+      - uses: actions/checkout@v6
       - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.julia-version }}
@@ -62,7 +62,7 @@ jobs:
         julia-arch:
           - x64
     steps:
-      - uses: actions/checkout@v5
+      - uses: actions/checkout@v6
       - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.julia-version }}
@@ -75,7 +75,7 @@ jobs:
   docs:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v5
+      - uses: actions/checkout@v6
       - uses: julia-actions/setup-julia@v2
         with:
           version: 'nightly'
diff --git a/src/higherorderfns.jl b/src/higherorderfns.jl
index 009928c3..24f96596 100644
--- a/src/higherorderfns.jl
+++ b/src/higherorderfns.jl
@@ -121,7 +121,7 @@ const SpBroadcasted2{Style<:SPVM,Axes,F,Args<:Tuple{SparseVecOrMat,SparseVecOrMa
 @inline columns(A::AbstractCompressedVector) = 1
 @inline columns(A::AbstractSparseMatrixCSC) = axes(A,2)
 @inline colrange(A::AbstractCompressedVector, j) = 1:length(nonzeroinds(A))
-@inline colrange(A::AbstractSparseMatrixCSC, j) = nzrange(A, j)
+Base.@propagate_inbounds colrange(A::AbstractSparseMatrixCSC, j) = nzrange(A, j)
 @inline colstartind(A::AbstractCompressedVector, j) = one(indtype(A))
 @inline colboundind(A::AbstractCompressedVector, j) = convert(indtype(A), length(nonzeroinds(A)) + 1)
 @inline colstartind(A::AbstractSparseMatrixCSC, j) = getcolptr(A)[j]
diff --git a/src/linalg.jl b/src/linalg.jl
index ad83fc51..0bc9f7ce 100644
--- a/src/linalg.jl
+++ b/src/linalg.jl
@@ -69,7 +69,7 @@ Base.@constprop :aggressive function spdensemul!(C, tA, tB, A, B, alpha, beta)
         T = eltype(C)
         _mul!(rangefun, diagop, odiagop, C, A, B, T(alpha), T(beta))
     else
-        @stable_muladdmul LinearAlgebra._generic_matmatmul!(C, 'N', 'N', wrap(A, tA), wrap(B, tB), MulAddMul(alpha, beta))
+        @stable_muladdmul LinearAlgebra._generic_matmatmul!(C, wrap(A, tA), wrap(B, tB), MulAddMul(alpha, beta))
     end
     return C
 end
@@ -232,7 +232,7 @@ function (*)(Da::Diagonal, A::SparseMatrixCSC, Db::Diagonal)
     rows = rowvals(A)
     vals = nonzeros(A)
     da, db = map(parent, (Da, Db))
-    for col in axes(A,2)
+    @inbounds for col in axes(A,2)
         dbcol = db[col]
         for i in nzrange(A, col)
             row = rows[i]
@@ -1351,7 +1351,7 @@ function _dot(x::AbstractSparseVector, A::AbstractSparseMatrixCSC, y::AbstractSp
         end
     end
     # diagonal
-    for i in axes(A,1)
+    @inbounds for i in axes(A,1)
         r1 = Int(Acolptr[i])
         r2 = Int(Acolptr[i+1]-1)
         r1 > r2 && continue
@@ -1391,7 +1391,7 @@ function ldiv!(D::Diagonal, A::Union{AbstractSparseMatrixCSC, AbstractSparseVect
     nonz = nonzeros(A)
     Arowval = rowvals(A)
     b = D.diag
-    for i=axes(b,1)
+    @inbounds for i=axes(b,1)
         iszero(b[i]) && throw(SingularException(i))
     end
     @inbounds for col in axes(A,2), p in nzrange(A, col)
@@ -1406,10 +1406,10 @@ function triu(S::AbstractSparseMatrixCSC{Tv,Ti}, k::Integer=0) where {Tv,Ti}
     m,n = size(S)
     colptr = Vector{Ti}(undef, n+1)
     nnz = 0
-    for col = 1 : min(max(k+1,1), n+1)
+    @inbounds for col = 1 : min(max(k+1,1), n+1)
         colptr[col] = 1
     end
-    for col = max(k+1,1) : n
+    @inbounds for col = max(k+1,1) : n
         for c1 in nzrange(S, col)
             rowvals(S)[c1] > col - k && break
             nnz += 1
@@ -1418,7 +1418,7 @@ function triu(S::AbstractSparseMatrixCSC{Tv,Ti}, k::Integer=0) where {Tv,Ti}
     end
     rowval = Vector{Ti}(undef, nnz)
     nzval = Vector{Tv}(undef, nnz)
-    for col = max(k+1,1) : n
+    @inbounds for col = max(k+1,1) : n
         c1 = getcolptr(S)[col]
         for c2 in colptr[col]:colptr[col+1]-1
             rowval[c2] = rowvals(S)[c1]
@@ -1434,7 +1434,7 @@ function tril(S::AbstractSparseMatrixCSC{Tv,Ti}, k::Integer=0) where {Tv,Ti}
     colptr = Vector{Ti}(undef, n+1)
     nnz = 0
     colptr[1] = 1
-    for col = 1 : min(n, m+k)
+    @inbounds for col = 1 : min(n, m+k)
         l1 = getcolptr(S)[col+1]-1
         for c1 = 0 : (l1 - getcolptr(S)[col])
             rowvals(S)[l1 - c1] < col - k && break
@@ -1442,12 +1442,12 @@ function tril(S::AbstractSparseMatrixCSC{Tv,Ti}, k::Integer=0) where {Tv,Ti}
         end
         colptr[col+1] = nnz+1
     end
-    for col = max(min(n, m+k)+2,1) : n+1
+    @inbounds for col = max(min(n, m+k)+2,1) : n+1
         colptr[col] = nnz+1
     end
     rowval = Vector{Ti}(undef, nnz)
     nzval = Vector{Tv}(undef, nnz)
-    for col = 1 : min(n, m+k)
+    @inbounds for col = 1 : min(n, m+k)
         c1 = getcolptr(S)[col+1]-1
         l2 = colptr[col+1]-1
         for c2 = 0 : l2 - colptr[col]
@@ -1469,8 +1469,8 @@ function sparse_diff1(S::AbstractSparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
     rowval = Vector{Ti}(undef, numnz)
     nzval = Vector{Tv}(undef, numnz)
     numnz = 0
-    colptr[1] = 1
-    for col = 1 : n
+    @inbounds colptr[1] = 1
+    @inbounds for col = 1 : n
         last_row = 0
         last_val = 0
         for k in nzrange(S, col)
@@ -1514,20 +1514,22 @@ function sparse_diff2(a::AbstractSparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
     rowval_a = rowvals(a)
     nzval_a = nonzeros(a)
 
-    ptrS = 1
-    colptr[1] = 1
+    @inbounds begin
+        ptrS = 1
+        colptr[1] = 1
 
-    n == 0 && return SparseMatrixCSC(m, n, colptr, rowval, nzval)
+        n == 0 && return SparseMatrixCSC(m, n, colptr, rowval, nzval)
 
-    startA = colptr_a[1]
-    stopA = colptr_a[2]
+        startA = colptr_a[1]
+        stopA = colptr_a[2]
 
-    rA = startA : stopA - 1
-    rowvalA = rowval_a[rA]
-    nzvalA = nzval_a[rA]
-    lA = stopA - startA
+        rA = startA : stopA - 1
+        rowvalA = rowval_a[rA]
+        nzvalA = nzval_a[rA]
+        lA = stopA - startA
+    end
 
-    for col = 1:n-1
+    @inbounds for col = 1:n-1
         startB, stopB = startA, stopA
         startA = colptr_a[col+1]
         stopA = colptr_a[col+2]
@@ -1614,7 +1616,7 @@ function opnorm(A::AbstractSparseMatrixCSC, p::Real=2)
         Tsum = promote_type(Float64,Tnorm)
         if p==1
             nA::Tsum = 0
-            for j in axes(A,2)
+            @inbounds for j in axes(A,2)
                 colSum::Tsum = 0
                 for i in nzrange(A, j)
                     colSum += abs(nonzeros(A)[i])
@@ -1626,7 +1628,7 @@ function opnorm(A::AbstractSparseMatrixCSC, p::Real=2)
             throw(ArgumentError("2-norm not yet implemented for sparse matrices. Try opnorm(Array(A)) or opnorm(A, p) where p=1 or Inf."))
         elseif p==Inf
             rowSum = zeros(Tsum,m)
-            for i in axes(nonzeros(A),1)
+            @inbounds for i in axes(nonzeros(A),1)
                 rowSum[rowvals(A)[i]] += abs(nonzeros(A)[i])
             end
             return convert(Tnorm, maximum(rowSum))
@@ -1683,8 +1685,8 @@ function opnormestinv(A::AbstractSparseMatrixCSC{T}, t::Integer = min(2,maximum(
 
     # Generate the block matrix
     X = Matrix{Ti}(undef, n, t)
-    X[1:n,1] .= 1
-    for j = 2:t
+    @inbounds X[1:n,1] .= 1
+    @inbounds for j = 2:t
         while true
             rand!(view(X,1:n,j), (-1, 1))
             yaux = X[1:n,j]' * X[1:n,1:j-1]
@@ -2019,8 +2021,8 @@ function mergeinds!(C::AbstractSparseMatrixCSC, A::AbstractSparseMatrixCSC)
     C_colptr = getcolptr(C)
     for col in axes(A,2)
         n_extra = 0
-        for ind in nzrange(A, col)
-            row = rowvals(A)[ind]
+        for ind in @inbounds nzrange(A, col)
+            row = @inbounds rowvals(A)[ind]
             row_exists, ind = rowcheck_index(C, row, col)
             if !row_exists
                 n_extra += 1
@@ -2052,31 +2054,31 @@ function mul!(C::AbstractSparseMatrixCSC, A::AbstractSparseMatrixCSC, D::Diagona
     if beta_is_zero || identical_nzinds
         identical_nzinds || copyinds!(C, A, copy_rows = !rows_match, copy_cols = !cols_match)
         resize!(Cnzval, length(Anzval))
-        if beta_is_zero
+        @inbounds if beta_is_zero
             if isone(alpha)
                 for col in axes(A,2), p in nzrange(A, col)
-                    @inbounds Cnzval[p] = Anzval[p] * b[col]
+                    Cnzval[p] = Anzval[p] * b[col]
                 end
             else
                 for col in axes(A,2), p in nzrange(A, col)
-                    @inbounds Cnzval[p] = Anzval[p] * b[col] * alpha
+                    Cnzval[p] = Anzval[p] * b[col] * alpha
                 end
             end
         else
             if isone(alpha)
                 for col in axes(A,2), p in nzrange(A, col)
-                    @inbounds Cnzval[p] = Anzval[p] * b[col] + Cnzval[p] * beta
+                    Cnzval[p] = Anzval[p] * b[col] + Cnzval[p] * beta
                 end
             else
                 for col in axes(A,2), p in nzrange(A, col)
-                    @inbounds Cnzval[p] = Anzval[p] * b[col] * alpha + Cnzval[p] * beta
+                    Cnzval[p] = Anzval[p] * b[col] * alpha + Cnzval[p] * beta
                 end
             end
         end
     else
         mergeinds!(C, A)
-        for col in axes(C,2), p in nzrange(C, col)
-            row = rowvals(C)[p]
+        for col in axes(C,2), p in @inbounds nzrange(C, col)
+            row = @inbounds rowvals(C)[p]
             # check if the index (row, col) is stored in A
             row_exists, row_ind_A = rowcheck_index(A, row, col)
             if row_exists
diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl
index 38a6748b..ae7292c7 100644
--- a/src/sparsematrix.jl
+++ b/src/sparsematrix.jl
@@ -217,13 +217,13 @@ julia> nnz(A)
 3
 ```
 """
-nnz(S::AbstractSparseMatrixCSC) = Int(getcolptr(S)[size(S, 2) + 1]) - 1
+nnz(S::AbstractSparseMatrixCSC) = @inbounds Int(getcolptr(S)[size(S, 2) + 1]) - 1
 nnz(S::ReshapedArray{<:Any,1,<:AbstractSparseMatrixCSC}) = nnz(parent(S))
 nnz(S::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}) = nnz(parent(S))
 nnz(S::UpperTriangular{<:Any,<:AbstractSparseMatrixCSC}) = nnz1(S)
 nnz(S::LowerTriangular{<:Any,<:AbstractSparseMatrixCSC}) = nnz1(S)
 nnz(S::SparseMatrixCSCColumnSubset) = nnz1(S)
-nnz1(S) = sum(length.(nzrange.(Ref(S), axes(S, 2))))
+nnz1(S) = @inbounds sum(length.(nzrange.(Ref(S), axes(S, 2))))
 
 function Base._simple_count(pred, S::AbstractSparseMatrixCSC, init::T) where T
     init + T(count(pred, nzvalview(S)) + pred(zero(eltype(S)))*(prod(size(S)) - nnz(S)))
@@ -308,8 +308,8 @@ of sparse array `A`. In conjunction with [`nonzeros`](@ref) and
 !!! warning
     Adding or removing nonzero elements to the matrix may invalidate the `nzrange`, one should not mutate the matrix while iterating.
 """
-nzrange(S::AbstractSparseMatrixCSC, col::Integer) = getcolptr(S)[col]:(getcolptr(S)[col+1]-1)
-nzrange(S::SparseMatrixCSCColumnSubset, col::Integer) = nzrange(S.parent, S.indices[2][col])
+Base.@propagate_inbounds nzrange(S::AbstractSparseMatrixCSC, col::Integer) = getcolptr(S)[col]:(getcolptr(S)[col+1]-1)
+Base.@propagate_inbounds nzrange(S::SparseMatrixCSCColumnSubset, col::Integer) = nzrange(S.parent, S.indices[2][col])
 nzrange(S::UpperTriangular{<:Any,<:SparseMatrixCSCUnion}, i::Integer) = nzrangeup(S.data, i)
 nzrange(S::LowerTriangular{<:Any,<:SparseMatrixCSCUnion}, i::Integer) = nzrangelo(S.data, i)
 
@@ -317,7 +317,7 @@ const AbstractSparseMatrixCSCInclAdjointAndTranspose = Union{AbstractSparseMatri
 function Base.isstored(A::AbstractSparseMatrixCSC, i::Integer, j::Integer)
     @boundscheck checkbounds(A, i, j)
     rows = rowvals(A)
-    for istored in nzrange(A, j) # could do binary search if the row indices are sorted?
+    @inbounds for istored in nzrange(A, j) # could do binary search if the row indices are sorted?
         i == rows[istored] && return true
     end
     return false
@@ -326,7 +326,7 @@ end
 function Base.isstored(A::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}, i::Integer, j::Integer)
     @boundscheck checkbounds(A, i, j)
     cols = rowvals(parent(A))
-    for istored in nzrange(parent(A), i)
+    @inbounds for istored in nzrange(parent(A), i)
         j == cols[istored] && return true
     end
     return false
@@ -674,7 +674,7 @@ function Base.copyto!(A::Array{T}, S::SparseMatrixCSC{<:Number}) where {T<:Numbe
     rowval = getrowval(S)
     nzval = getnzval(S)
     linear_index_col0 = 0   # Linear index before column (linear index = linear_index_col0 + row)
-    for col in axes(S, 2)
+    @inbounds for col in axes(S, 2)
         for i in nzrange(S, col)
             row = rowval[i]
             val = nzval[i]
@@ -2269,7 +2269,7 @@ function (+)(A::SparseMatrixCSCUnion, B::Array)
     C = Ref(zero(eltype(A))) .+ B
     rowinds, nzvals = rowvals(A), nonzeros(A)
     for j in axes(A,2)
-        for i in nzrange(A, j)
+        @inbounds for i in nzrange(A, j)
             rowidx = rowinds[i]
             C[rowidx,j] = nzvals[i] + B[rowidx,j]
         end
@@ -2281,7 +2281,7 @@ function (+)(A::Array, B::SparseMatrixCSCUnion)
     C = A .+ Ref(zero(eltype(B)))
     rowinds, nzvals = rowvals(B), nonzeros(B)
     for j in axes(B,2)
-        for i in nzrange(B, j)
+        @inbounds for i in nzrange(B, j)
             rowidx = rowinds[i]
             C[rowidx,j] = A[rowidx,j] + nzvals[i]
         end
@@ -2293,7 +2293,7 @@ function (-)(A::SparseMatrixCSCUnion, B::Array)
     C = Ref(zero(eltype(A))) .- B
     rowinds, nzvals = rowvals(A), nonzeros(A)
     for j in axes(A,2)
-        for i in nzrange(A, j)
+        @inbounds for i in nzrange(A, j)
             rowidx = rowinds[i]
             C[rowidx,j] = nzvals[i] - B[rowidx,j]
         end
@@ -2305,7 +2305,7 @@ function (-)(A::Array, B::SparseMatrixCSCUnion)
     C = A .- Ref(zero(eltype(B)))
     rowinds, nzvals = rowvals(B), nonzeros(B)
     for j in axes(B,2)
-        for i in nzrange(B, j)
+        @inbounds for i in nzrange(B, j)
             rowidx = rowinds[i]
             C[rowidx,j] = A[rowidx,j] - nzvals[i]
         end
@@ -4086,7 +4086,7 @@ function is_hermsym(A::AbstractSparseMatrixCSC, check::Function)
     rowval = rowvals(A)
     nzval = nonzeros(A)
     tracker = copy(getcolptr(A))
-    for col in axes(A,2)
+    @inbounds for col in axes(A,2)
         # `tracker` is updated such that, for symmetric matrices,
         # the loop below starts from an element at or below the
         # diagonal element of column `col`"
@@ -4165,7 +4165,7 @@ function istriu(A::AbstractSparseMatrixCSC, k::Integer=0)
     rowval = rowvals(A)
     nzval  = nonzeros(A)
 
-    for col = 1:min(n, m-1)
+    @inbounds for col = 1:min(n, m-1)
         l1 = colptr[col+1]-1
         for i = 0 : (l1 - colptr[col])
             if rowval[l1-i] <= col - k
@@ -4186,7 +4186,7 @@ function istril(A::AbstractSparseMatrixCSC, k::Integer=0)
     rowval = rowvals(A)
     nzval  = nonzeros(A)
 
-    for col = 2:n
+    @inbounds for col = 2:n
         for i = colptr[col] : (colptr[col+1]-1)
             if rowval[i] >= col - k
                 # subsequent rows would also lie below the band
@@ -4603,7 +4603,7 @@ function copytrito!(M::AbstractMatrix, S::AbstractSparseMatrixCSC, uplo::Char)
 
     rv = rowvals(S)
     nz = nonzeros(S)
-    for col in axes(S,2)
+    @inbounds for col in axes(S,2)
         trirange = uplo == 'U' ? (1:min(col, size(S,1))) : (col:size(S,1))
         fill!(view(M, trirange, col), zero(eltype(S)))
         for i in nzrange(S, col)
diff --git a/src/sparsevector.jl b/src/sparsevector.jl
index 6734b986..6e47de67 100644
--- a/src/sparsevector.jl
+++ b/src/sparsevector.jl
@@ -115,7 +115,7 @@ end
 nnz(x::AbstractCompressedVector) = length(nonzeros(x))
 function nnz(x::SparseColumnView)
     rowidx, colidx = parentindices(x)
-    return length(nzrange(parent(x), colidx))
+    return length(@inbounds nzrange(parent(x), colidx))
 end
 nnz(x::SparseVectorView) = nnz(x.parent)
 nnz(x::SparseVectorPartialView) = length(nonzeroinds(x))
diff --git a/test/linalg.jl b/test/linalg.jl
index ce4f4e31..1e849305 100644
--- a/test/linalg.jl
+++ b/test/linalg.jl
@@ -109,6 +109,14 @@ end
     end
 end
 
+@testset "multiplication of special sparse with dense matrix" begin
+    # this results in a call of the most generic multiplication code in LinearAlgebra.jl
+    A = randn(2, 2)
+    S = sparse(A)
+    B = rand(1, 2)'
+    @test Symmetric(S) * B ≈ Symmetric(A) * B
+end
+
 @testset "sparse transpose adjoint" begin
     A = sprand(10, 10, 0.75)
     @test A' == SparseMatrixCSC(A')
diff --git a/test/runtests.jl b/test/runtests.jl
index c515e78d..5819db1a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,13 +7,22 @@ if Base.get_bool_env("SPARSEARRAYS_AQUA_TEST", false)
     include("ambiguous.jl")
 end
 
-for file in readlines(joinpath(@__DIR__, "testgroups"))
-    file == "" && continue # skip empty lines
-    include(file * ".jl")
-end
+include("allowscalar.jl")
+include("fixed.jl")
+include("higherorderfns.jl")
+include("sparsematrix_constructors_indexing.jl")
+include("sparsematrix_ops.jl")
+include("sparsevector.jl")
+include("issues.jl")
 
 if Base.USE_GPL_LIBS
 
+    include("cholmod.jl")
+    include("umfpack.jl")
+    include("spqr.jl")
+    include("linalg.jl")
+    include("linalg_solvers.jl")
+
     nt = @static if isdefined(Threads, :maxthreadid)
         Threads.maxthreadid()
     else
diff --git a/test/testgroups b/test/testgroups
deleted file mode 100644
index a547762c..00000000
--- a/test/testgroups
+++ /dev/null
@@ -1,12 +0,0 @@
-allowscalar
-cholmod
-fixed
-higherorderfns
-issues
-linalg
-linalg_solvers
-sparsematrix_constructors_indexing
-sparsematrix_ops
-sparsevector
-spqr
-umfpack