Merge branch 'main' into dk/cleanup-matprod_dest

Author: dkarrasch

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, 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

@@ -53,12 +53,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
@@ -108,7 +108,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)
@@ -119,6 +119,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) ||
@@ -1733,7 +1736,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/cholmod.jl

@@ -1287,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,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)))
 
@@ -154,7 +153,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) =
@@ -518,8 +516,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/sparsematrix.jl

@@ -313,6 +313,8 @@ Base.@propagate_inbounds nzrange(S::SparseMatrixCSCColumnSubset, col::Integer) =
 nzrange(S::UpperTriangular{<:Any,<:SparseMatrixCSCUnion}, i::Integer) = nzrangeup(S.data, i)
 nzrange(S::LowerTriangular{<:Any,<:SparseMatrixCSCUnion}, i::Integer) = nzrangelo(S.data, i)
 
+indtype(S::SparseMatrixCSCColumnSubset{<:Any,Ti}) where {Ti} = Ti
+
 const AbstractSparseMatrixCSCInclAdjointAndTranspose = Union{AbstractSparseMatrixCSC,Adjoint{<:Any,<:AbstractSparseMatrixCSC},Transpose{<:Any,<:AbstractSparseMatrixCSC}}
 function Base.isstored(A::AbstractSparseMatrixCSC, i::Integer, j::Integer)
     @boundscheck checkbounds(A, i, j)

test/linalg.jl

@@ -236,6 +236,20 @@ end
     end
 end
 
+@testset "Dense times symmetric/Hermitian sparse matrix multiplication" begin
+    A = [1 3; 2 4]
+    As = sparse(A)
+    B = [1 1; 1 1]
+    @test mul!(copy(B), B, Hermitian(A), true, true) == mul!(copy(B), B, Hermitian(As), true, true)
+end
+
+@testset "Column view of sparse matrix " begin
+    S = sparse(1:4, 1:4, 1:4)
+    Sv = @view S[:,3:4]
+    @test Sv * sparse(ones(2)) == Sv*ones(2) == Matrix(Sv) * ones(2)
+    @test Sv * sparse(ones(2,2)) == Sv*ones(2,2) == Matrix(Sv) * ones(2,2)
+end
+
 @testset "in-place sparse-sparse mul!" begin
     for n in (20, 30)
         sA = sprandn(ComplexF64, n, n, 0.1); A = Array(sA)

test/linalg_solvers.jl

@@ -130,20 +130,24 @@ end
 @testset "sparse matrix cond" begin
     Random.seed!(1235)
     local A = sparse(reshape([1.0], 1, 1))
-    Ac = sprandn(20, 20,.5) + im*sprandn(20, 20,.5)
-    Ar = sprandn(20, 20,.5) + eps()*I
     @test cond(A, 1) == 1.0
-    # For a discussion of the tolerance, see #14778
-    @test 0.99 <= cond(Ar, 1) \ opnorm(Ar, 1) * opnorm(inv(Array(Ar)), 1) < 3
-    @test 0.99 <= cond(Ac, 1) \ opnorm(Ac, 1) * opnorm(inv(Array(Ac)), 1) < 3
-    @test 0.99 <= cond(Ar, Inf) \ opnorm(Ar, Inf) * opnorm(inv(Array(Ar)), Inf) < 3
-    @test 0.99 <= cond(Ac, Inf) \ opnorm(Ac, Inf) * opnorm(inv(Array(Ac)), Inf) < 3
     @test_throws ArgumentError cond(A,2)
     @test_throws ArgumentError cond(A,3)
     Arect = spzeros(10, 6)
     @test_throws DimensionMismatch cond(Arect, 1)
     @test_throws ArgumentError cond(Arect,2)
     @test_throws DimensionMismatch cond(Arect, Inf)
+    Ac = sprandn(20, 20,.5) + im*sprandn(20, 20,.5)
+    Ar = sprandn(20, 20,.5) + eps()*I
+    # For a discussion of the tolerance, see #14778
+    @test 0.99 <= cond(Ar, 1) \ opnorm(Ar, 1) * opnorm(inv(Array(Ar)), 1) < 3
+    @test 0.99 <= cond(Ac, 1) \ opnorm(Ac, 1) * opnorm(inv(Array(Ac)), 1) < 3
+    @test 0.99 <= cond(Ar, Inf) \ opnorm(Ar, Inf) * opnorm(inv(Array(Ar)), Inf) < 3
+    @test 0.99 <= cond(Ac, Inf) \ opnorm(Ac, Inf) * opnorm(inv(Array(Ac)), Inf) < 3
+    #issue 680
+    A22 = sparse(randn(2,2))
+    @test 0.99 ≤ cond(Array(A22), 1) / cond(A22, 1) < 3
+    @test 0.99 ≤ cond(Array(A22), Inf) / cond(A22, Inf) < 3
 end
 
 @testset "sparse matrix opnormestinv" begin
@@ -159,6 +163,9 @@ end
     @test_throws ArgumentError SparseArrays.opnormestinv(Ac,0)
     @test_throws ArgumentError SparseArrays.opnormestinv(Ac,21)
     @test_throws DimensionMismatch SparseArrays.opnormestinv(sprand(3,5,.9))
+    #issue 680
+    A33 = sparse(randn(3,3))
+    @test SparseArrays.opnormestinv(A33,3) ≈ opnorm(inv(Array(A33)),1) atol=1e-4
 end
 
 @testset "factorization" begin