diff --git a/Project.toml b/Project.toml index dc33978..8ad7940 100644 --- a/Project.toml +++ b/Project.toml @@ -10,9 +10,9 @@ OhMyThreads = "67456a42-1dca-4109-a031-0a68de7e3ad5" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [compat] -julia = "1" LinearMaps = "3" OhMyThreads = "0.8.3" +julia = "1" [extras] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" diff --git a/README.md b/README.md index f66e456..b109690 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,24 @@ + + # BlockSparseMatrices +

+ Description +

+ [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://djukic14.github.io/BlockSparseMatrices.jl/stable/) [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://djukic14.github.io/BlockSparseMatrices.jl/dev/) [![Build Status](https://github.com/djukic14/BlockSparseMatrices.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/djukic14/BlockSparseMatrices.jl/actions/workflows/CI.yml?query=branch%3Amain) [![Coverage](https://codecov.io/gh/djukic14/BlockSparseMatrices.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/djukic14/BlockSparseMatrices.jl) + + +## Introduction +This package provides an AbstractMatrix object that allows to construct sparse block matrices. + +## Installation +The package can be installed by + +```@julia +import Pkg +Pkg.add("https://github.com/djukic14/BlockSparseMatrices.jl") +``` \ No newline at end of file diff --git a/docs/make.jl b/docs/make.jl index 03e80cf..7f91a7c 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -14,7 +14,18 @@ makedocs(; edit_link="main", assets=String[], ), - pages=["Home" => "index.md"], + pages=[ + "Introduction" => "index.md", + "Manual" => "manual.md", + "Further Details" => "details.md", + "Contributing" => "contributing.md", + "API Reference" => "apiref.md", + ], ) -deploydocs(; repo="github.com/djukic14/BlockSparseMatrices.jl", devbranch="main") +deploydocs(; + repo="github.com/djukic14/BlockSparseMatrices.jl", + target="build", + #devbranch="dev", + versions=["stable" => "v^"],#=, "dev" => "dev"=# +) diff --git a/docs/src/apiref.md b/docs/src/apiref.md new file mode 100644 index 0000000..e5add5b --- /dev/null +++ b/docs/src/apiref.md @@ -0,0 +1,14 @@ +# API Reference + +## Types +```@docs +BlockSparseMatrices.BlockSparseMatrix +BlockSparseMatrices.SymmetricBlockMatrix +BlockSparseMatrices.DenseMatrixBlock +``` + +## Functions +```@docs +BlockSparseMatrices.findcolor! +BlockSparseMatrices.islessinordering +``` \ No newline at end of file diff --git a/docs/src/assets/logo.svg b/docs/src/assets/logo.svg new file mode 100644 index 0000000..4599475 --- /dev/null +++ b/docs/src/assets/logo.svg @@ -0,0 +1,191 @@ + +image/svg+xml \ No newline at end of file diff --git a/docs/src/assets/nocolor.svg b/docs/src/assets/nocolor.svg new file mode 100644 index 0000000..8d9cff1 --- /dev/null +++ b/docs/src/assets/nocolor.svg @@ -0,0 +1,123 @@ + +image/svg+xml \ No newline at end of file diff --git a/docs/src/contributing.md b/docs/src/contributing.md new file mode 100644 index 0000000..9189893 --- /dev/null +++ b/docs/src/contributing.md @@ -0,0 +1,55 @@ +# Contributing + +In order to contribute to this package directly create a pull request against the `main` branch. Before doing so please: + +- Follow the style of the surrounding code. +- Supplement the documentation. +- Write tests and check that no errors occur. + + +--- +## Style + +For a consistent style the [JuliaFormatter.jl](https://github.com/domluna/JuliaFormatter.jl) package is used which enforces the style defined in the *.JuliaFormatter.toml* file. To follow this style simply run +```julia +using JuliaFormatter +format(pkgdir(BlockSparseMatrices)) +``` + +!!! note + That all files follow the JuliaFormatter style is tested during the unit tests. Hence, do not forget to execute the two lines above. Otherwise, the tests are likely to not pass. + + +--- +## Documentation + +Add documentation for any changes or new features following the style of the existing documentation. For more information you can have a look at the [Documenter.jl](https://documenter.juliadocs.org/stable/) documentation. + + +--- +## [Tests](@id tests) + +Write tests for your code changes and verify that no errors occur, e.g., by running +```julia +using Pkg +Pkg.test("BlockSparseMatrices.jl") +``` +For a detailed information on which parts are tested the coverage can be evaluated on your local machine, e.g., by +```julia +using Pkg +Pkg.test("BlockSparseMatrices"; coverage=true, julia_args=["-t 4"]) + +# determine coverage +using Coverage +src_folder = pkgdir(BlockSparseMatrices) * "/src" +coverage = process_folder(src_folder) +LCOV.writefile("path-to-folder-you-like" * "BlockSparseMatrices.lcov.info", coverage) + +clean_folder(src_folder) # delete .cov files + +# extract information about coverage +covered_lines, total_lines = get_summary(coverage) +@info "Current coverage:\n$covered_lines of $total_lines lines ($(round(Int, covered_lines / total_lines * 100)) %)" +``` + +In Visual Studio Code the [Coverage Gutters](https://marketplace.visualstudio.com/items?itemName=ryanluker.vscode-coverage-gutters) plugin can be used to visualize the tested lines of the code by inserting the path of the *BlockSparseMatrices.lcov.info* file in the settings. \ No newline at end of file diff --git a/docs/src/details.md b/docs/src/details.md new file mode 100644 index 0000000..be64559 --- /dev/null +++ b/docs/src/details.md @@ -0,0 +1,36 @@ + +# Further Details + +## Matrix-Vector Product +For the multithreaded matrix-vector product + +$$\bm y = \bm A \bm x,$$ + +where $\bm A$ is a `BlockSparseMatrix` or a `SymmetricBlockMatrix`, $\bm x$ a random vector, and $\bm y$ the solution vector, we use a coloring scheme to allow a multithreaded matrix-vector product without unnecessary allocations. +Each block is assigned to one color such that all blocks of the same color do not share common row or column indices. +This allows to evaluate all matrix-vector products within one color parallel without allocating multiple solution vectors $\bm y$. + +The coloring scheme works best for large numbers of dense blocks. + +The row and column indices of blocks should not be overlapping. Only unique sets of row and column indices are supported as well as row and column indices that are subsets of the indices of other blocks. + +### Example +Example of the coloring scheme using four groups. Each color can be evaluated using multiple threads with a single solution vector + +```@raw html + + +
+ +
+ +
Figure 1: Dense blocks in a block sparse matrix.
+
+ +
+ +
Figure 2: Colored dense blocks.
+
+ +
+``` diff --git a/docs/src/index.md b/docs/src/index.md index 1c09b22..d4fcbbf 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,14 +1,15 @@ -```@meta -CurrentModule = BlockSparseMatrices -``` +# BlockSparseMatrices.jl -# BlockSparseMatrices +This package provides an abstract matrix type for block structured matrices, where the majority of the matrix is zero. +This is in general the case for the near interactions in fast methods for the boundary element method. -Documentation for [BlockSparseMatrices](https://github.com/djukic14/BlockSparseMatrices.jl). +--- +## Installation -```@index -``` +Installing BlockSparseMatrices is done by entering the package manager (enter `]` at the julia REPL) and issuing: -```@autodocs -Modules = [BlockSparseMatrices] ``` +pkg> add BlockSparseMatrices +``` + + diff --git a/docs/src/manual.md b/docs/src/manual.md new file mode 100644 index 0000000..df14480 --- /dev/null +++ b/docs/src/manual.md @@ -0,0 +1,30 @@ + +# General Usage + +The basic usage of this package is provided in the following examples. + +## Example: Block Sparse Matrix +```@example blocksparsematrix +using BlockSparseMatrices + +# define dense matrix blocks +mat1 = randn(ComplexF64, 2, 2) +block1 = BlockSparseMatrices.DenseMatrixBlock(mat1, 1:2, 1:2) +mat2 = randn(ComplexF64, 3, 3) +block2 = BlockSparseMatrices.DenseMatrixBlock(mat2, 3:5, 3:5) + +blockmatrix = BlockSparseMatrix([block1, block2], 5, 5) +``` + +## Example: Symmetric Block Sparse Matrix +```@example symblocksparsematrix +using BlockSparseMatrices + +# define dense matrix blocks +mat1 = randn(ComplexF64, 2, 2) +diagonalblock = BlockSparseMatrices.DenseMatrixBlock(mat1, 1:2, 1:2) +mat2 = randn(ComplexF64, 3, 3) +offdiagonalblock = BlockSparseMatrices.DenseMatrixBlock(mat2, 7:9, 7:9) + +diagonalmatrix = SymmetricBlockMatrix([diagonalblock], [offdiagonalblock], 9, 9) +``` \ No newline at end of file diff --git a/src/blockmatrix.jl b/src/blockmatrix.jl index 1b6c269..013cbb7 100644 --- a/src/blockmatrix.jl +++ b/src/blockmatrix.jl @@ -1,3 +1,19 @@ +""" + struct BlockSparseMatrix{T,M,D} <: AbstractBlockMatrix{T} + +Sparse block matrix. + +# Fields +- `blocks::Vector{M}`: Dense blocks. +- `size::Tuple{Int,Int}`: Size of the matrix. +- `forwardbuffer::Vector{T}`: ???. +- `adjointbuffer::Vector{T}`: ???. +- `buffer::Vector{T}`: ???. +- `rowindexdict::D`: Global rowindices of dense blocks. +- `colindexdict::D`: Global columnindices of dense blocks. +- `threadsafecolors::Vector{Vector{Int}}`: Coloring of dense blocks for multithreaded matrix-vector product. +- `ntasks::Int`: Threads used for the matrix Vector product. +""" struct BlockSparseMatrix{T,M,D} <: AbstractBlockMatrix{T} blocks::Vector{M} size::Tuple{Int,Int} diff --git a/src/matrixblock/abstractmatrixblock.jl b/src/matrixblock/abstractmatrixblock.jl index 9231057..800047a 100644 --- a/src/matrixblock/abstractmatrixblock.jl +++ b/src/matrixblock/abstractmatrixblock.jl @@ -88,6 +88,11 @@ function Base.axes(block::AbstractMatrixBlock) end #TODO: check if this is the correct ordering or if it should be reversed -> benchmark performance +""" + islessinordering(blocka::AbstractMatrixBlock, blockb::AbstractMatrixBlock) + +Sorting rule for `AbstractMatrixBlock` objects +""" function islessinordering(blocka::AbstractMatrixBlock, blockb::AbstractMatrixBlock) if maximum(rowindices(blocka)) < maximum(rowindices(blockb)) return true @@ -100,6 +105,7 @@ struct isthreadsafe end struct issymthreadsafe isthreadsafe::isthreadsafe end + issymthreadsafe() = issymthreadsafe(isthreadsafe()) function (ists::issymthreadsafe)(blocka, blockb) @@ -128,6 +134,24 @@ function (::isthreadsafe)(blocka, blockb) end end +""" + findcolor!( + blockid::Int, + threadsafecolors::AbstractArray, + blocks::Vector{M}; + threadsafecheck=isthreadsafe(), + color=1, + ) where {M} + +Assigns recursively threadsafe color to the dense block with index `blockid`. + +# Arguments +- `blockid::Int`: Index of block. +- `threadsafecolors::AbstractArray`: Contains assigned blocks. +- `blocks::Vector{M}`: Vector with all blocks. +- `threadsafecheck=isthreadsafe()`: Threadsafe check, differs between symmetric and non-symmetric matrix. +- `color=1`: Currently tested color, increased if block does not fit in color. +""" function findcolor!( blockid::Int, threadsafecolors::AbstractArray, diff --git a/src/matrixblock/densematrixblock.jl b/src/matrixblock/densematrixblock.jl index eba313c..02f7b31 100644 --- a/src/matrixblock/densematrixblock.jl +++ b/src/matrixblock/densematrixblock.jl @@ -1,3 +1,13 @@ +""" + DenseMatrixBlock{T,M,RC} <: AbstractMatrixBlock{T} + +Dense block in a sparse blockmatrix. + +# Fields +- `matrix::M`: Matrix. +- `rowindices::RC`: Global row indices. +- `colindices::RC`: Global column indices. +""" struct DenseMatrixBlock{T,M,RC} <: AbstractMatrixBlock{T} matrix::M rowindices::RC diff --git a/src/symmetricblockmatrix.jl b/src/symmetricblockmatrix.jl index 2fd10a6..2464be0 100644 --- a/src/symmetricblockmatrix.jl +++ b/src/symmetricblockmatrix.jl @@ -1,3 +1,22 @@ +""" + SymmetricBlockMatrix{T,DM,M,D} <: AbstractBlockMatrix{T} + +Symmetric sparse block matrix. + +# Fields +- `diagonals::Vector{DM}`: Dense blocks on the diagonal of the matrix. +- `offdiagonals::Vector{M}`: Dense blocks not on the diagonal. +- `size::Tuple{Int,Int}`: Size of the matrix. +- `forwardbuffer::Vector{T}`: ???. +- `adjointbuffer::Vector{T}`: ???. +- `buffer::Vector{T}`: ???. +- `diagonalrowindexdict::D`: Global rowindices of diagonal dense blocks. +- `diagonalcolindexdict::D`: Global columnindices of offdiagonal dense blocks. +- `offdiagonalrowindexdict::D`: Global rowindices of diagonal dense blocks. +- `offdiagonalcolindexdict::D`: Global columnindices of offdiagonal dense blocks. +- `threadsafecolors::Vector{Vector{Int}}`: Coloring of dense blocks for multithreaded matrix-vector product. +- `ntasks::Int`: Threads used for the matrix Vector product. +""" struct SymmetricBlockMatrix{T,DM,M,D} <: AbstractBlockMatrix{T} diagonals::Vector{DM} offdiagonals::Vector{M} diff --git a/test/test_blockmatrix.jl b/test/test_blockmatrix.jl index 06b5dca..ac29422 100644 --- a/test/test_blockmatrix.jl +++ b/test/test_blockmatrix.jl @@ -26,11 +26,14 @@ block5 = BlockSparseMatrices.DenseMatrixBlock(mat1, 6:7, 6:7) @test !BlockSparseMatrices.isthreadsafe()(block4, block1) blockmatrix = BlockSparseMatrix([block1, block2, block3, block4], 5, 5) +@test blockmatrix[1, 1] == mat1[1, 1] blockmatrix2 = BlockSparseMatrix([block1, block2, block3, block4], (5, 5)) blockmatrix3 = BlockSparseMatrix( [mat1, mat2, mat3, mat4], [1:2, 3:5, 1:2, 3:5], [1:2, 3:5, 3:5, 1:2], (5, 5) ) blockmatrix4 = BlockSparseMatrix([block1, block2, block5], 7, 7) +@test blockmatrix4[1, 5] == 0.0 + @test length(blockmatrix4.threadsafecolors) == 2 @test blockmatrix.blocks == blockmatrix2.blocks == blockmatrix3.blocks @@ -63,3 +66,25 @@ LinearAlgebra.mul!(x2, adjoint(M), x, α, β) LinearAlgebra.mul!(x1, transpose(blockmatrix), x, α, β) LinearAlgebra.mul!(x2, transpose(M), x, α, β) @test x1 ≈ x2 + +blockmatrix[1, 1] = 1.0 +@test blockmatrix[1, 1] == 1.0 +@test BlockSparseMatrices.eachblockindex(blockmatrix) == + BlockSparseMatrices.eachblockindex(transpose(blockmatrix)) + +block1 = BlockSparseMatrices.DenseMatrixBlock(mat1, 1:2, 1:2) +block2 = BlockSparseMatrices.DenseMatrixBlock(mat2, 3:5, 3:5) +block3 = BlockSparseMatrices.DenseMatrixBlock(mat3, 1:2, 3:5) +block4 = BlockSparseMatrices.DenseMatrixBlock(mat4, 3:5, 1:2) +block5 = BlockSparseMatrices.DenseMatrixBlock(mat1, 6:7, 6:7) + +@test BlockSparseMatrices.isthreadsafe()(block1, block2) +@test BlockSparseMatrices.isthreadsafe()(block2, block1) +@test BlockSparseMatrices.isthreadsafe()(block3, block4) +@test BlockSparseMatrices.isthreadsafe()(block4, block3) +@test !BlockSparseMatrices.isthreadsafe()(block1, block3) +@test !BlockSparseMatrices.isthreadsafe()(block3, block1) +@test !BlockSparseMatrices.isthreadsafe()(block1, block4) +@test !BlockSparseMatrices.isthreadsafe()(block4, block1) + +blockmatrix = BlockSparseMatrix([block1, block2, block3, block4], 5, 5) diff --git a/test/test_symmetricblockmatrix.jl b/test/test_symmetricblockmatrix.jl index d8d103e..27a0ca3 100644 --- a/test/test_symmetricblockmatrix.jl +++ b/test/test_symmetricblockmatrix.jl @@ -30,6 +30,12 @@ diagonalmatrix3 = SymmetricBlockMatrix( @test size(adjoint(diagonalmatrix)) == (5, 5) @test size(adjoint(diagonalmatrix), 1) == size(adjoint(diagonalmatrix), 2) == 5 @test nnz(diagonalmatrix) == 25 +@test diagonalmatrix[1, 1] == block1[1, 1] +@test diagonalmatrix[3, 3] == block2[1, 1] +@test diagonalmatrix[1, 3] == block3[1, 1] +@test diagonalmatrix[3, 1] == block3[1, 1] +@test BlockSparseMatrices.eachoffdiagonalindex(diagonalmatrix) == + BlockSparseMatrices.eachoffdiagonalindex(transpose(diagonalmatrix)) x = rand(ComplexF64, 5) @test diagonalmatrix * x ≈ M * x @@ -52,6 +58,13 @@ LinearAlgebra.mul!(x1, transpose(diagonalmatrix), x, α, β) LinearAlgebra.mul!(x2, transpose(M), x, α, β) @test x1 ≈ x2 +diagonalmatrix[1, 1] = 1+im*2 +@test diagonalmatrix[1, 1] == 1+im*2 +diagonalmatrix[1, 3] = 1+im*2 +@test diagonalmatrix[1, 3] == 1+im*2 +diagonalmatrix[3, 1] = 2+im*1 +@test diagonalmatrix[3, 1] == 2+im*1 + # check threadsafecolors block1 = randn(ComplexF64, 2, 2) block2 = randn(ComplexF64, 2, 2)