From 9e20eaf64efee93278aa5417d10c573a92e86d91 Mon Sep 17 00:00:00 2001 From: Philippe Maincon Date: Tue, 15 Aug 2023 09:13:45 +0200 Subject: [PATCH] blocks can be added separately to a big sparse matrix --- src/Muscade.jl | 2 +- src/SparseCat.jl | 107 +++++++++++++++++++++++++++++++++--------- test/TestSparseCat.jl | 23 ++++++--- 3 files changed, 102 insertions(+), 30 deletions(-) diff --git a/src/Muscade.jl b/src/Muscade.jl index 6b7a26a..cb60fb9 100644 --- a/src/Muscade.jl +++ b/src/Muscade.jl @@ -51,7 +51,7 @@ module Muscade export solve include("SparseCat.jl") - export cat! + export blocksparse,cat!,addin!,zero! include("StaticX.jl") diff --git a/src/SparseCat.jl b/src/SparseCat.jl index 3b055b3..039b4d7 100644 --- a/src/SparseCat.jl +++ b/src/SparseCat.jl @@ -16,66 +16,94 @@ # irow = ilr[ilv] """ - bigsparse = cat(blocks::Matrix{<:SparseMatrixCSC}) + bigsparse,start = blocksparse(blocks::Matrix{<:SparseMatrixCSC}) -Assemble sparse blocs into a large sparse matrix. +Prepare for the assembly of sparse blocs into a large sparse matrix. Unassigned blocks (undef'ed) will be treated as all-zero. +`bigsparse` is allocated, with the correct sparsity structure, but its `nzval` undef'ed. +Where some blocks share the same sparsity structure, `blocks` can have `===` elements. -See also: [`cat!`](@ref) +See also: [`zero!`](@ref),[`addin!`](@ref),[`cat!`](@ref) """ -function Base.cat(B::Matrix{SparseMatrixCSC{Tv,Ti}}) where{Tv,Ti<:Integer} +function blocksparse(B::Matrix{SparseMatrixCSC{Tv,Ti}}) where{Tv,Ti<:Integer} nbr,nbc = size(B) nbr>0 && nbc>0 || muscadeerror("must have length(B)>0") pgr = 𝕫1(undef,nbr+1) # pointers into solution vector, with global*solution=rhs pgr[1] = 1 for ibr = 1:nbr - pgr[ibr+1] = pgr[ibr]+size(B[ibr,1],1) + nlr = 0 + for ibc = 1:nbc + if isassigned(B,ibr,ibc) + nlr = size(B[ibr,ibc],1) + break + end + end + nlr > 0 || muscadeerror("every bloc-row must contain at least one assigned Sparse") + pgr[ibr+1] = pgr[ibr]+nlr end ngr = pgr[end]-1 pgc = 𝕫1(undef,nbc+1) # pointers into rhs vector pgc[1] = 1 for ibc = 1:nbc - pgc[ibc+1] = pgc[ibc]+size(B[1,ibc],2) + nlc = 0 + for ibr = 1:nbr + if isassigned(B,ibr,ibc) + nlc = size(B[ibr,ibc],2) + break + end + end + nlc > 0 || muscadeerror("every bloc-column must contain at least one assigned Sparse") + pgc[ibc+1] = pgc[ibc]+nlc end ngc = pgc[end]-1 ngv = 0 for ibc = 1:nbc for ibr = 1:nbr - ngv += nnz(B[ibr,ibc]) + if isassigned(B,ibr,ibc) + ngv += nnz(B[ibr,ibc]) + end end end - igr = 𝕫1(undef,ngv) - gv = 𝕣1(undef,ngv) + igr = 𝕫1(undef,ngv ) + gv = 𝕣1(undef,ngv ) pigr = 𝕫1(undef,ngc+1) pigr[1] = 1 + start = Vector{𝕫2}(undef,nbr) + for ibc = 1:nbc + nlc = pgc[ibc+1]-pgc[ibc] + start[ibc] = zeros(𝕫,nbr,nlc) + end + igv = 1 for ibc = 1:nbc for ilc = 1:size(B[1,ibc],2) igc = ilc + pgc[ibc]-1 for ibr = 1:nbr - pilr,ilr,lv = B[ibr,ibc].colptr, B[ibr,ibc].rowval,B[ibr,ibc].nzval - for ilv = pilr[ilc]:pilr[ilc+1]-1 - igr[igv]= pgr[ibr]-1 + ilr[ilv] - gv[igv] = lv[ilv] - igv += 1 + if isassigned(B,ibr,ibc) + start[ibc][ibr,ilc] = igv + pilr,ilr = B[ibr,ibc].colptr, B[ibr,ibc].rowval + for ilv = pilr[ilc]:pilr[ilc+1]-1 + igr[igv]= pgr[ibr]-1 + ilr[ilv] + igv += 1 + end end pigr[igc+1] = igv end end end - return SparseMatrixCSC(ngr,ngc,pigr,igr,gv) + return SparseMatrixCSC(ngr,ngc,pigr,igr,gv),start end """ - bigsparse = cat!(bigsparse,blocs::Matrix{<:SparseMatrixCSC}) + cat!(bigsparse,blocs::Matrix{<:SparseMatrixCSC}) -Assemble sparse blocs into a large sparse matrix. Will fail silently or throw an error unless +Assemble sparse blocs into a preallocated sparse matrix. Will fail silently or throw an error unless `bigsparse` has the correct sparsity structure for the given `blocs`. -See also: [`cat`](@ref) +See also: [`blocksparse!`](@ref),[`zero!`](@ref),[`addin!`](@ref) """ function cat!(out::SparseMatrixCSC{Tv,Ti},B::Matrix{SparseMatrixCSC{Tv,Ti}}) where{Tv,Ti<:Integer} nbr,nbc = size(B) @@ -83,12 +111,47 @@ function cat!(out::SparseMatrixCSC{Tv,Ti},B::Matrix{SparseMatrixCSC{Tv,Ti}}) whe for ibc = 1:nbc for ilc = 1:size(B[1,ibc],2) for ibr = 1:nbr - pilr,lv = B[ibr,ibc].colptr,B[ibr,ibc].nzval - for ilv = pilr[ilc]:pilr[ilc+1]-1 - gv[igv] = lv[ilv] - igv += 1 + if isassigned(B,ibr,ibc) + pilr,lv = B[ibr,ibc].colptr,B[ibr,ibc].nzval + for ilv = pilr[ilc]:pilr[ilc+1]-1 + gv[igv] = lv[ilv] + igv += 1 + end end end end end end +""" + addin!(bigsparse,blocs::Matrix{<:SparseMatrixCSC},start,ibr,ibc) + +Add a sparse into one of the blocks of a large sparse matrix. Will fail silently or throw an error unless +`bigsparse` has the correct sparsity structure for the given `blocs`. Use [`sparseblocks`](@ref) to + create `bigsparse` and `start`. + +See also: [`blocksparse!`](@ref),[`zero!`](@ref),[`cat!`](@ref) +""" +function addin!(out::SparseMatrixCSC{Tv,Ti},B::SparseMatrixCSC{Tv,Ti},start,ibr::𝕫,ibc::𝕫) where{Tv,Ti<:Integer} + gv = out.nzval + start[ibc][ibr,1] > 0 || muscadeerror("Trying to cat! into an empty block") + for ilc = 1:size(B,2) + igv = start[ibc][ibr,ilc] + pilr,lv = B.colptr,B.nzval + for ilv = pilr[ilc]:pilr[ilc+1]-1 + gv[igv]+= lv[ilv] + igv += 1 + end + end +end +""" + zero!(sparse) + +Set to zero the vector `nzval` of values in a sparse matrix + +See also: [`blocksparse!`](@ref),[`addin!`](@ref),[`cat!`](@ref) +""" +function zero!(out::SparseMatrixCSC) + out.nzval .= 0 +end + + diff --git a/test/TestSparseCat.jl b/test/TestSparseCat.jl index edc6508..fe6b90b 100644 --- a/test/TestSparseCat.jl +++ b/test/TestSparseCat.jl @@ -10,19 +10,28 @@ ncol = 2 B = Matrix{SparseMatrixCSC{𝕣,𝕫}}(undef,nrow,ncol) for irow = 1:nrow for icol = 1:ncol - B[irow,icol] = sparse([1,1,2,3,3],[1,2,2,2,3],randn(5)) + if irow