Merge branch 'main' into export_abstractmatrixcsc

Author: j-fu

src/linalg.jl

@@ -205,6 +205,27 @@ const SparseOrTri{Tv,Ti} = Union{SparseMatrixCSCUnion{Tv,Ti},SparseTriangular{Tv
 *(A::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}, B::SparseOrTri) = spmatmul(copy(A), B)
 *(A::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}, B::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}) = spmatmul(copy(A), copy(B))
 
+(*)(Da::Diagonal, A::Union{SparseMatrixCSCUnion, AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}}, Db::Diagonal) = Da * (A * Db)
+function (*)(Da::Diagonal, A::SparseMatrixCSC, Db::Diagonal)
+    (size(Da, 2) == size(A,1) && size(A,2) == size(Db,1)) ||
+        throw(DimensionMismatch("incompatible sizes"))
+    T = promote_op(matprod, eltype(Da), promote_op(matprod, eltype(A), eltype(Db)))
+    dest = similar(A, T)
+    vals_dest = nonzeros(dest)
+    rows = rowvals(A)
+    vals = nonzeros(A)
+    da, db = map(parent, (Da, Db))
+    for col in axes(A,2)
+        dbcol = db[col]
+        for i in nzrange(A, col)
+            row = rows[i]
+            val = vals[i]
+            vals_dest[i] = da[row] * val * dbcol
+        end
+    end
+    dest
+end
+
 # Gustavson's matrix multiplication algorithm revisited.
 # The result rowval vector is already sorted by construction.
 # The auxiliary Vector{Ti} xb is replaced by a Vector{Bool} of same length.
@@ -629,6 +650,62 @@ function dot(A::AbstractSparseMatrixCSC, B::Union{DenseMatrixUnion,WrapperMatrix
     return conj(dot(B, A))
 end
 
+function dot(x::AbstractSparseVector, D::Diagonal, y::AbstractVector)
+    d = D.diag
+    if length(x) != length(y) || length(y) != length(d)
+        throw(
+            DimensionMismatch("Vectors and matrix have different dimensions, x has a length $(length(x)), y has a length $(length(y)), D has side dimension $(length(d))")
+        )
+    end
+    nzvals = nonzeros(x)
+    nzinds = nonzeroinds(x)
+    s = zero(typeof(dot(first(x), first(D), first(y))))
+    @inbounds for nzidx in eachindex(nzvals)
+        s += dot(nzvals[nzidx], d[nzinds[nzidx]], y[nzinds[nzidx]])
+    end
+    return s
+end
+
+dot(x::AbstractVector, D::Diagonal, y::AbstractSparseVector) = adjoint(dot(y, D', x))
+
+function dot(x::AbstractSparseVector, D::Diagonal, y::AbstractSparseVector)
+    d = D.diag
+    if length(y) != length(x) || length(y) != length(d)
+        throw(
+            DimensionMismatch("Vectors and matrix have different dimensions, x has a length $(length(x)), y has a length $(length(y)), Q has side dimension $(length(d))")
+        )
+    end
+    xnzind = nonzeroinds(x)
+    ynzind = nonzeroinds(y)
+    xnzval = nonzeros(x)
+    ynzval = nonzeros(y)
+    s = zero(typeof(dot(first(x), first(D), first(y))))
+    if isempty(xnzind) || isempty(ynzind)
+        return s
+    end
+
+    x_idx = 1
+    y_idx = 1
+    x_idx_last = length(xnzind)
+    y_idx_last = length(ynzind)
+
+    # go through the nonzero indices of a and b simultaneously
+    @inbounds while x_idx <= x_idx_last && y_idx <= y_idx_last
+        ix = xnzind[x_idx]
+        iy = ynzind[y_idx]
+        if ix == iy
+            s += dot(xnzval[x_idx], d[ix], ynzval[y_idx])
+            x_idx += 1
+            y_idx += 1
+        elseif ix < iy
+            x_idx += 1
+        else
+            y_idx += 1
+        end
+    end
+    return s
+end
+
 ## triangular sparse handling
 ## triangular multiplication
 function LinearAlgebra.generic_trimatmul!(C::StridedVecOrMat, uploc, isunitc, tfun::Function, A::SparseMatrixCSCUnion, B::AbstractVecOrMat)

src/solvers/spqr.jl

@@ -146,7 +146,7 @@ Matrix{T}(Q::QRSparseQ) where {T} = lmul!(Q, Matrix{T}(I, size(Q, 1), min(size(Q
 
 # From SPQR manual p. 6
 _default_tol(A::AbstractSparseMatrixCSC) =
-    20*sum(size(A))*eps(real(eltype(A)))*maximum(norm(view(A, :, i)) for i in 1:size(A, 2))
+    20*sum(size(A))*eps()*maximum(norm(view(A, :, i)) for i in 1:size(A, 2))
 
 """
     qr(A::SparseMatrixCSC; tol=_default_tol(A), ordering=ORDERING_DEFAULT) -> QRSparse

test/linalg.jl

@@ -850,6 +850,20 @@ end
             @test dot(TA,WB) ≈ dot(Matrix(TA), WB)
             @test dot(TA,WC) ≈ dot(Matrix(TA), WC)
         end
+        for M in (A, B, C)
+            D = Diagonal(M * M')
+            a = spzeros(Complex{Float64}, size(D, 1))
+            a[1:3] = rand(Complex{Float64}, 3)
+            b = spzeros(Complex{Float64}, size(D, 1))
+            b[1:3] = rand(Complex{Float64}, 3)
+            @test dot(a, D, b) ≈ dot(a, sparse(D), b)
+            @test dot(b, D, a) ≈ dot(b, sparse(D), a)
+            @test dot(b, D, a) ≈ dot(b, D, collect(a))
+            @test dot(b, D, a) ≈ dot(collect(b), D, a)
+            @test_throws DimensionMismatch dot(b, D, [a; 1])
+            @test_throws DimensionMismatch dot([b; 1], D, a)
+            @test_throws DimensionMismatch dot([b; 1], D, [a; 1])
+        end
     end
     @test_throws DimensionMismatch dot(sprand(5,5,0.2),sprand(5,6,0.2))
     @test_throws DimensionMismatch dot(rand(5,5),sprand(5,6,0.2))
@@ -912,4 +926,24 @@ end
     @test sparse(3I, 4, 5) == sparse(1:4, 1:4, 3, 4, 5)
     @test sparse(3I, 5, 4) == sparse(1:4, 1:4, 3, 5, 4)
 end
+
+@testset "diagonal-sandwiched triple multiplication" begin
+    S = sprand(4, 6, 0.2)
+    D1 = Diagonal(axes(S,1))
+    D2 = Diagonal(axes(S,2) .+ 4)
+    A = Array(S)
+    C = D1 * S * D2
+    @test C isa SparseMatrixCSC
+    @test C ≈ D1 * A * D2
+    C = D2 * S' * D1
+    @test C isa SparseMatrixCSC
+    @test C ≈ D2 * A' * D1
+    C = D1 * view(S, :, :) * D2
+    @test C isa SparseMatrixCSC
+    @test C ≈ D1 * A * D2
+
+    @test_throws DimensionMismatch D2 * S * D2
+    @test_throws DimensionMismatch D1 * S * D1
+end
+
 end