From afe05b5d4e329bd4b78cf3a545250bbaec9ec387 Mon Sep 17 00:00:00 2001 From: Petr Krysl Date: Thu, 19 Mar 2026 11:22:06 -0700 Subject: [PATCH] add transpose and adjoint --- NEWS.md | 4 ++ Project.toml | 2 +- src/SparseMatricesCSR.jl | 2 +- src/SparseMatrixCSR.jl | 83 ++++++++++++++++++++++++++++++++++++---- test/SparseMatrixCSR.jl | 20 ++++++++++ 5 files changed, 102 insertions(+), 9 deletions(-) diff --git a/NEWS.md b/NEWS.md index e2c67cf..209688f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,6 +5,10 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [0.6.13] - 2026.03.19 + +- Transpose and adjoint added. + ## [0.6.12] - 2026.03.12 - Minor fixes in documentation. diff --git a/Project.toml b/Project.toml index d25aa2c..43625f1 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SparseMatricesCSR" uuid = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1" authors = ["Víctor Sande ", "Francesc Verdugo "] -version = "0.6.12" +version = "0.6.13" [deps] Atomix = "a9b6321e-bd34-4604-b9c9-b65b8de01458" diff --git a/src/SparseMatricesCSR.jl b/src/SparseMatricesCSR.jl index 4ea4ef3..7f4b9d3 100644 --- a/src/SparseMatricesCSR.jl +++ b/src/SparseMatricesCSR.jl @@ -4,7 +4,7 @@ using SparseArrays using LinearAlgebra using SuiteSparse -import Base: convert, copy, size, getindex, setindex!, show, count, *, IndexStyle +import Base: convert, copy, size, getindex, setindex!, show, count, *, IndexStyle, transpose import LinearAlgebra: mul!, lu, lu! import SparseArrays: nnz, getnzval, nonzeros, nzvalview, nzrange import SparseArrays: findnz, rowvals, getnzval, issparse diff --git a/src/SparseMatrixCSR.jl b/src/SparseMatrixCSR.jl index 838edbf..9e038a3 100644 --- a/src/SparseMatrixCSR.jl +++ b/src/SparseMatrixCSR.jl @@ -134,6 +134,73 @@ function Base.copy(a::SparseMatrixCSR{Bi}) where Bi SparseMatrixCSR{Bi}(a.m,a.n,copy(a.rowptr),copy(a.colval),copy(a.nzval)) end +widelength(x::SparseMatrixCSR) = prod(Int64.(size(x))) + +function Base.sizehint!(S::SparseMatrixCSR, n::Integer) + nhint = min(n, widelength(S)) + sizehint!(getcolval(S), nhint) + sizehint!(nonzeros(S), nhint) + return S +end + +Base.adjoint(A::SparseMatrixCSR) = Adjoint(A) +Base.transpose(A::SparseMatrixCSR) = Transpose(A) +Base.copy(A::Transpose{<:Any,<:SparseMatrixCSR}) = + ftranspose_csr(A.parent, identity) +Base.copy(A::Adjoint{<:Any,<:SparseMatrixCSR}) = + ftranspose_csr(A.parent, adjoint) + +function ftranspose_csr(A::SparseMatrixCSR{Bi,TvA,Ti}, f::Function, eltype::Type{Tv} = TvA) where {Bi,Tv,TvA,Ti} + X = SparseMatrixCSR{Bi}(size(A, 2), size(A, 1), + ones(Ti, size(A, 2)+1), + Vector{Ti}(undef, 0), + Vector{Tv}(undef, 0)) + sizehint!(X, nnz(A)) + return halfperm_csr!(X, A, axes(A,1), f) +end + +function halfperm_csr!(X::SparseMatrixCSR{Bi,Tv,Ti}, A::SparseMatrixCSR{Bi,TvA,Ti}, + q::AbstractVector{<:Integer}, f::F = identity) where {Bi,Tv,TvA,Ti,F<:Function} + _computecolptrs_halfperm_csr!(X, A) + _distributevals_halfperm_csr!(X, A, q, f) + return X +end + +function _computecolptrs_halfperm_csr!(X::SparseMatrixCSR{Bi,Tv,Ti}, A::SparseMatrixCSR{Bi,TvA,Ti}) where {Bi,Tv,TvA,Ti} + # Compute transpose(A[:,q])'s row counts. Store shifted forward one position in X.rowptr. + fill!(X.rowptr, 0) + o = getoffset(A) + @inbounds for k in 1:nnz(A) + X.rowptr[A.colval[k] + o + 1] += 1 + end + # Compute transpose(A[:,q])'s row pointers. Store shifted forward one position in X.rowptr. + X.rowptr[1] = Bi + countsum = Bi + @inbounds for k in 2:(size(A, 2) + 1) + overwritten = X.rowptr[k] + X.rowptr[k] = countsum + countsum += overwritten + end +end + +function _distributevals_halfperm_csr!(X::SparseMatrixCSR{Bi,Tv,Ti}, + A::SparseMatrixCSR{Bi,TvA,Ti}, q::AbstractVector{<:Integer}, f::F) where {Bi,Tv,TvA,Ti,F<:Function} + resize!(X.nzval, nnz(A)) + resize!(X.colval, nnz(A)) + o = getoffset(A) + @inbounds for Xi in axes(A,1) + Aj = q[Xi] + for Ak in nzrange(A, Aj) + Ai = A.colval[Ak] + o + Xk = X.rowptr[Ai + 1] + X.colval[Xk] = Xi - 1 + Bi # since colval are Bi-based + X.nzval[Xk] = f(A.nzval[Ak]) + X.rowptr[Ai + 1] += 1 + end + end + return +end + _copy_and_increment(x) = copy(x) .+ 1 function LinearAlgebra.fillstored!(a::SparseMatrixCSR,v) @@ -220,6 +287,8 @@ getrowptr(S::SparseMatrixCSR) = S.rowptr getnzval(S::SparseMatrixCSR) = S.nzval getcolval(S::SparseMatrixCSR) = S.colval +colvals(S::SparseMatrixCSR) = S.colval + """ getBi(S::SparseMatrixCSR{Bi}) where {Bi} @@ -263,7 +332,7 @@ nonzeros(S::SparseMatrixCSR) = S.nzval nzvalview(S::SparseMatrixCSR) = view(nonzeros(S), 1:nnz(S)) """ - colvals(S::SparseMatrixCSR{Bi}) where {Bi} + colvals(S::SparseMatrixCSR) Return a vector of the col indices of `S`. The stored values are indexes to arrays with `Bi`-based indexing, but the `colvals(S)` array itself is a standard 1-based @@ -273,15 +342,15 @@ Providing access to how the col indices are stored internally can be useful in conjunction with iterating over structural nonzero values. See also [`nonzeros`](@ref) and [`nzrange`](@ref). """ -colvals(S::SparseMatrixCSR) = S.colval -""" - nzrange(S::SparseMatrixCSR{Bi}, row::Integer) where {Bi} -Return the range of indices to the structural nonzero values of a -sparse matrix row. The returned range of indices is always 1-based even for `Bi != 1`. """ -nzrange(S::SparseMatrixCSR{Bi}, row::Integer) where {Bi} = S.rowptr[row]+getoffset(S):S.rowptr[row+1]-Bi + nzrange(S::SparseMatrixCSR{Bi}, row::Integer) + + Return the range of indices to the structural nonzero values of a + sparse matrix row. The returned range of indices is always 1-based even for `Bi != 1`. +""" + nzrange(S::SparseMatrixCSR{Bi}, row::Integer) where {Bi} = S.rowptr[row]+getoffset(S):S.rowptr[row+1]-Bi """ findnz(S::SparseMatrixCSR{Bi,Tv,Ti}) diff --git a/test/SparseMatrixCSR.jl b/test/SparseMatrixCSR.jl index 35b7403..ee2b366 100644 --- a/test/SparseMatrixCSR.jl +++ b/test/SparseMatrixCSR.jl @@ -139,6 +139,26 @@ function test_csr(Bi,Tv,Ti) @test nnz(A) == length(SparseArrays.nzvalview(A)) == 3 @test SparseArrays.nzvalview(A) == [4., 5, 6] end + + # Test transpose + CSC = sparse(I,J,V,maxrows,maxcols) + if Bi == 1 + CSR = sparsecsr(I,J,V,maxrows,maxcols) + @test CSR' == CSC' + @test copy(CSR') == CSC' + @test CSR' == transpose(CSR) + @test copy(CSR') == transpose(CSR) + @test CSR' == adjoint(CSR) + @test typeof(copy(CSR')) == typeof(copy(transpose(CSR))) == SparseMatrixCSR{Bi,Tv,Ti} + end + # Rectangular matrix transpose + A = Tv[7 0 -3 0 -1 0; + 2 8 0 0 0 0; + 0 0 1 0 0 0; + 0 0 0 -2 0 6] + As = sparsecsr(A) + @test As' == transpose(As) == adjoint(As) == sparse(A') + @test typeof(copy(As')) == typeof(copy(transpose(As))) end function test_lu(Bi,I,J,V)