diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml new file mode 100644 index 00000000..a1a04cb0 --- /dev/null +++ b/.buildkite/pipeline.yml @@ -0,0 +1,18 @@ +env: + SECRET_CODECOV_TOKEN: "l7PxpMgxHEiz2l3UDDj5122Xirv5du1QFQiTf09YISC8VOz/nXeQC2yYY6iTMbccQG3xrumfWke9vtLbAc9LPH75kw3zQODuWMbz0eGSf9wtw7foqO6KAo5neMTeZJ0ZFgFMa2Q89T5SNIQsi4b5Zs7BVto6qB9Z/3QEs/BpsR25cYkY4Y6JBU6XuqZ6GRAc6BtlB4OmBEp3BBsXatx6y64zF0qbp4rocmcBQEoeQkcXtxf0dfA0KNAHpweWPrzIAMZ0aUYp6iEikxwLY5TjVFhOIcpVXUlxIdhl2qaaDF6b6WlgGXiGLpBjsLQfvlOgIlXH59Ddg0IVxboF2a37OA==;U2FsdGVkX19MveWMe1aasYZoJwaeiKO5XqwMZ/utOFj7h1CE9UX4nduBMXTyJ77tpNmCBYsQNJ/PadbX2224hQ==" + +steps: + - label: "Julia v1" + plugins: + - JuliaCI/julia#v1: + version: "1" + - JuliaCI/julia-test#v1: ~ + - JuliaCI/julia-coverage#v1: + codecov: true + agents: + queue: "juliagpu" + cuda: "*" + if: build.message !~ /\[skip tests\]/ + timeout_in_minutes: 60 + env: + JULIA_SMC_TEST_GROUP: 'GPU' \ No newline at end of file diff --git a/.github/workflows/Benchmark.yml b/.github/workflows/Benchmark.yml index fda51ed6..4a4a444b 100644 --- a/.github/workflows/Benchmark.yml +++ b/.github/workflows/Benchmark.yml @@ -12,11 +12,11 @@ jobs: runs-on: ubuntu-latest if: contains(github.event.pull_request.labels.*.name, 'benchmark') steps: - - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v2 + - uses: actions/checkout@v6 + - uses: julia-actions/setup-julia@v3 with: version: "1" - - uses: julia-actions/cache@v2 + - uses: julia-actions/cache@v3 - name: Extract Package Name from Project.toml id: extract-package-name run: | @@ -46,14 +46,14 @@ jobs: echo '' >> body.md echo '' >> body.md - name: Find Comment - uses: peter-evans/find-comment@v3 + uses: peter-evans/find-comment@v4 id: fcbenchmark with: issue-number: ${{ github.event.pull_request.number }} comment-author: 'github-actions[bot]' body-includes: Benchmark Results - name: Comment on PR - uses: peter-evans/create-or-update-comment@v4 + uses: peter-evans/create-or-update-comment@v5 with: comment-id: ${{ steps.fcbenchmark.outputs.comment-id }} issue-number: ${{ github.event.pull_request.number }} diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml index 1f74cd19..492a6c0b 100644 --- a/.github/workflows/CompatHelper.yml +++ b/.github/workflows/CompatHelper.yml @@ -15,7 +15,7 @@ jobs: run: which julia continue-on-error: true - name: Install Julia, but only if it is not already available in the PATH - uses: julia-actions/setup-julia@v2 + uses: julia-actions/setup-julia@v3 with: version: '1' arch: ${{ runner.arch }} diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml index a8d1110e..2f2e57df 100644 --- a/.github/workflows/Documentation.yml +++ b/.github/workflows/Documentation.yml @@ -19,11 +19,11 @@ jobs: statuses: write runs-on: ubuntu-latest steps: - - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v2 + - uses: actions/checkout@v6 + - uses: julia-actions/setup-julia@v3 with: version: '1' - - uses: julia-actions/cache@v2 + - uses: julia-actions/cache@v3 - name: Install dependencies run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' - name: Build and deploy diff --git a/.github/workflows/Format.yml b/.github/workflows/Format.yml new file mode 100644 index 00000000..dd2b2c50 --- /dev/null +++ b/.github/workflows/Format.yml @@ -0,0 +1,12 @@ +name: Format suggestions +on: + pull_request: + # this argument is not required if you don't use the `suggestion-label` input + types: [opened, reopened, synchronize, labeled, unlabeled] +jobs: + code-style: + runs-on: ubuntu-latest + steps: + - uses: julia-actions/julia-format@v4 + with: + version: '1' # Set `version` to '1.0.54' if you need to use JuliaFormatter.jl v1.0.54 (default: '1') diff --git a/.github/workflows/Test-GPU.yml b/.github/workflows/Test-GPU.yml deleted file mode 100644 index cd67c6e1..00000000 --- a/.github/workflows/Test-GPU.yml +++ /dev/null @@ -1,46 +0,0 @@ -name: Test-GPU - -on: - push: - branches: - - main - pull_request: - -concurrency: - group: ${{ github.workflow }}-${{ github.ref }} - cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }} - -# needed to allow julia-actions/cache to delete old caches that it has created -permissions: - actions: write - contents: read - -jobs: - test: - runs-on: self-hosted - env: - CUDA_VISIBLE_DEVICES: 1 - JULIA_DEPOT_PATH: /scratch/github-actions/julia_depot_alexis - JULIA_SMC_TEST_GROUP: "GPU" - strategy: - matrix: - julia-version: ['1.10', '1'] - steps: - - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v2 - with: - version: ${{ matrix.julia-version }} - arch: x64 - - uses: julia-actions/julia-downgrade-compat@v1 - if: ${{ matrix.version == '1.10' }} - with: - skip: LinearAlgebra, Random, SparseArrays - - uses: julia-actions/cache@v2 - - uses: julia-actions/julia-buildpkg@v1 - - uses: julia-actions/julia-runtest@v1 - - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v5 - with: - files: lcov.info - token: ${{ secrets.CODECOV_TOKEN }} - fail_ci_if_error: false diff --git a/.github/workflows/Test.yml b/.github/workflows/Test.yml index 1f0046a1..6e8c6af5 100644 --- a/.github/workflows/Test.yml +++ b/.github/workflows/Test.yml @@ -20,24 +20,24 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - julia-version: ['1.10', '1'] + julia-version: ['1.10', '1.11', '1.12'] steps: - - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v2 + - uses: actions/checkout@v6 + - uses: julia-actions/setup-julia@v3 with: version: ${{ matrix.julia-version }} arch: x64 - - uses: julia-actions/julia-downgrade-compat@v1 + - uses: julia-actions/julia-downgrade-compat@v2 if: ${{ matrix.version == '1.10' }} with: skip: LinearAlgebra, Random, SparseArrays - - uses: julia-actions/cache@v2 + - uses: julia-actions/cache@v3 - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v5 + - uses: codecov/codecov-action@v6 with: files: lcov.info token: ${{ secrets.CODECOV_TOKEN }} - fail_ci_if_error: false \ No newline at end of file + fail_ci_if_error: false diff --git a/CITATION.cff b/CITATION.cff index c2c0a8fb..a36bc54e 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -18,7 +18,7 @@ identifiers: - type: doi value: 10.5281/zenodo.11314275 description: Zenodo -repository-code: 'https://github.com/gdalle/SparseMatrixColorings.jl' +repository-code: 'https://github.com/JuliaDiff/SparseMatrixColorings.jl' abstract: >- Coloring algorithms for sparse Jacobian and Hessian matrices diff --git a/Project.toml b/Project.toml index 6656e125..72a5fad2 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SparseMatrixColorings" uuid = "0a514795-09f3-496d-8182-132a7b665d35" +version = "0.4.27" authors = ["Guillaume Dalle", "Alexis Montoison"] -version = "0.4.21" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" @@ -12,23 +12,41 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [weakdeps] +BandedMatrices = "aae01518-5342-5314-be14-df237901396f" +BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" +BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8" Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" +GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7" +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" +MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" +cuSPARSE = "b26da814-b3bc-49ef-b0ee-c816305aa060" [extensions] -SparseMatrixColoringsCUDAExt = "CUDA" +SparseMatrixColoringsBandedMatricesExt = "BandedMatrices" +SparseMatrixColoringsBlockBandedMatricesExt = ["BlockArrays", "BlockBandedMatrices"] +SparseMatrixColoringsCUDAExt = ["CUDA", "cuSPARSE"] SparseMatrixColoringsCliqueTreesExt = "CliqueTrees" SparseMatrixColoringsColorsExt = "Colors" +SparseMatrixColoringsGPUArraysExt = "GPUArrays" +SparseMatrixColoringsJuMPExt = ["JuMP", "MathOptInterface"] [compat] ADTypes = "1.2.1" -CUDA = "5.8.2" +BandedMatrices = "1.9.4" +BlockArrays = "1.6.3" +BlockBandedMatrices = "0.13.1" +CUDA = "6.0.0" CliqueTrees = "1" Colors = "0.12.11, 0.13" DocStringExtensions = "0.8,0.9" -LinearAlgebra = "<0.0.1, 1" +GPUArrays = "11.5.0" +JuMP = "1.29.1" +LinearAlgebra = "1" +MathOptInterface = "1.45.0" PrecompileTools = "1.2.1" -Random = "<0.0.1, 1" -SparseArrays = "<0.0.1, 1" +Random = "1" +SparseArrays = "1" +cuSPARSE = "6.0.0" julia = "1.10" diff --git a/README.md b/README.md index 66b10dd3..b37b0ca8 100644 --- a/README.md +++ b/README.md @@ -1,9 +1,11 @@ # SparseMatrixColorings.jl -[![Build Status](https://github.com/gdalle/SparseMatrixColorings.jl/actions/workflows/Test.yml/badge.svg?branch=main)](https://github.com/gdalle/SparseMatrixColorings.jl/actions/workflows/Test.yml?query=branch%3Amain) -[![Stable Documentation](https://img.shields.io/badge/docs-stable-blue.svg)](https://gdalle.github.io/SparseMatrixColorings.jl/stable/) -[![Dev Documentation](https://img.shields.io/badge/docs-dev-blue.svg)](https://gdalle.github.io/SparseMatrixColorings.jl/dev/) -[![Coverage](https://codecov.io/gh/gdalle/SparseMatrixColorings.jl/branch/main/graph/badge.svg)](https://app.codecov.io/gh/gdalle/SparseMatrixColorings.jl) +[![Build Status](https://github.com/juliadiff/SparseMatrixColorings.jl/actions/workflows/Test.yml/badge.svg?branch=main)](https://github.com/juliadiff/SparseMatrixColorings.jl/actions/workflows/Test.yml?query=branch%3Amain) +[![GPU build status](https://badge.buildkite.com/7d8ed289d7bdb5a25ae48b2c778a202ce4990b7ee558cdfef8.svg?branch=main)](https://buildkite.com/julialang/sparsematrixcolorings-dot-jl) +[![Coverage](https://codecov.io/gh/juliadiff/SparseMatrixColorings.jl/branch/main/graph/badge.svg)](https://app.codecov.io/gh/juliadiff/SparseMatrixColorings.jl) + +[![Stable Documentation](https://img.shields.io/badge/docs-stable-blue.svg)](https://juliadiff.org/SparseMatrixColorings.jl/stable/) +[![Dev Documentation](https://img.shields.io/badge/docs-dev-blue.svg)](https://juliadiff.org/SparseMatrixColorings.jl/dev/) [![Code Style: Blue](https://img.shields.io/badge/code%20style-blue-4495d1.svg)](https://github.com/JuliaDiff/BlueStyle) [![arXiv](https://img.shields.io/badge/arXiv-2505.07308-b31b1b.svg)](https://arxiv.org/abs/2505.07308) [![DOI](https://zenodo.org/badge/801999408.svg)](https://zenodo.org/doi/10.5281/zenodo.11314275) diff --git a/docs/make.jl b/docs/make.jl index 78fba268..bbf3f859 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -2,7 +2,10 @@ using Documenter using DocumenterInterLinks using SparseMatrixColorings -links = InterLinks("ADTypes" => "https://sciml.github.io/ADTypes.jl/stable/") +links = InterLinks( + "ADTypes" => "https://sciml.github.io/ADTypes.jl/stable/", + "BandedMatrices" => "https://julialinearalgebra.github.io/BandedMatrices.jl/stable/", +) cp(joinpath(@__DIR__, "..", "README.md"), joinpath(@__DIR__, "src", "index.md"); force=true) @@ -21,5 +24,7 @@ makedocs(; ) deploydocs(; - repo="github.com/gdalle/SparseMatrixColorings.jl", push_preview=true, devbranch="main" + repo="github.com/JuliaDiff/SparseMatrixColorings.jl", + push_preview=true, + devbranch="main", ) diff --git a/docs/src/api.md b/docs/src/api.md index c121d1e3..640198d3 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -17,8 +17,15 @@ SparseMatrixColorings coloring fast_coloring ColoringProblem +``` + +## Coloring algorithms + +```@docs GreedyColoringAlgorithm ConstantColoringAlgorithm +OptimalColoringAlgorithm +StructuredColoringAlgorithm ``` ## Result analysis diff --git a/docs/src/dev.md b/docs/src/dev.md index 2f253dbe..cd294a17 100644 --- a/docs/src/dev.md +++ b/docs/src/dev.md @@ -41,6 +41,7 @@ SparseMatrixColorings.TreeSetColoringResult SparseMatrixColorings.LinearSystemColoringResult SparseMatrixColorings.BicoloringResult SparseMatrixColorings.remap_colors +SparseMatrixColorings.decompress_csc! ``` ## Testing @@ -50,6 +51,9 @@ SparseMatrixColorings.directly_recoverable_columns SparseMatrixColorings.symmetrically_orthogonal_columns SparseMatrixColorings.structurally_orthogonal_columns SparseMatrixColorings.structurally_biorthogonal +SparseMatrixColorings.substitutable_columns +SparseMatrixColorings.substitutable_bidirectional +SparseMatrixColorings.rank_nonzeros_from_trees SparseMatrixColorings.valid_dynamic_order ``` @@ -58,7 +62,7 @@ SparseMatrixColorings.valid_dynamic_order ```@docs SparseMatrixColorings.respectful_similar SparseMatrixColorings.matrix_versions -SparseMatrixColorings.same_pattern +SparseMatrixColorings.compatible_pattern ``` ## Visualization @@ -76,3 +80,9 @@ SparseMatrixColorings.what_fig_61 SparseMatrixColorings.efficient_fig_1 SparseMatrixColorings.efficient_fig_4 ``` + +## Misc + +```@docs +SparseMatrixColorings.cycle_range +``` diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md index fec54f13..6e803511 100644 --- a/docs/src/tutorial.md +++ b/docs/src/tutorial.md @@ -31,7 +31,7 @@ problem = ColoringProblem() The algorithm defines how you want to solve it. It can be either a [`GreedyColoringAlgorithm`](@ref) or a [`ConstantColoringAlgorithm`](@ref). For `GreedyColoringAlgorithm`, you can select options such as -- the order in which vertices are processed (a subtype of [`AbstractOrder`](@ref SparseMatrixColorings.AbstractOrder)) +- the order in which vertices are processed (a subtype of [`AbstractOrder`](@ref SparseMatrixColorings.AbstractOrder) , or a tuple of such objects) - the type of decompression you want (`:direct` or `:substitution`) ```@example tutorial diff --git a/ext/SparseMatrixColoringsBandedMatricesExt.jl b/ext/SparseMatrixColoringsBandedMatricesExt.jl new file mode 100644 index 00000000..05d9677c --- /dev/null +++ b/ext/SparseMatrixColoringsBandedMatricesExt.jl @@ -0,0 +1,45 @@ +module SparseMatrixColoringsBandedMatricesExt + +using BandedMatrices: BandedMatrix, bandrange, bandwidths, colrange, rowrange +using SparseMatrixColorings: + BipartiteGraph, + ColoringProblem, + ColumnColoringResult, + StructuredColoringAlgorithm, + RowColoringResult, + column_colors, + cycle_range, + row_colors +import SparseMatrixColorings as SMC + +#= +This code is partly taken from ArrayInterface.jl and FiniteDiff.jl +https://github.com/JuliaArrays/ArrayInterface.jl +https://github.com/JuliaDiff/FiniteDiff.jl +=# + +function SMC.coloring( + A::BandedMatrix, + ::ColoringProblem{:nonsymmetric,:column}, + ::StructuredColoringAlgorithm; + kwargs..., +) + width = length(bandrange(A)) + color = cycle_range(width, size(A, 2)) + bg = BipartiteGraph(A) + return ColumnColoringResult(A, bg, color) +end + +function SMC.coloring( + A::BandedMatrix, + ::ColoringProblem{:nonsymmetric,:row}, + ::StructuredColoringAlgorithm; + kwargs..., +) + width = length(bandrange(A)) + color = cycle_range(width, size(A, 1)) + bg = BipartiteGraph(A) + return RowColoringResult(A, bg, color) +end + +end diff --git a/ext/SparseMatrixColoringsBlockBandedMatricesExt.jl b/ext/SparseMatrixColoringsBlockBandedMatricesExt.jl new file mode 100644 index 00000000..aaecb0ad --- /dev/null +++ b/ext/SparseMatrixColoringsBlockBandedMatricesExt.jl @@ -0,0 +1,97 @@ +module SparseMatrixColoringsBlockBandedMatricesExt + +using BlockArrays: blockaxes, blockfirsts, blocklasts, blocksize, blocklengths +using BlockBandedMatrices: + BandedBlockBandedMatrix, + BlockBandedMatrix, + blockbandrange, + blockbandwidths, + blocklengths, + blocksize, + subblockbandwidths +using SparseMatrixColorings: + BipartiteGraph, + ColoringProblem, + ColumnColoringResult, + StructuredColoringAlgorithm, + RowColoringResult, + column_colors, + cycle_range, + row_colors +import SparseMatrixColorings as SMC + +#= +This code is partly taken from ArrayInterface.jl and FiniteDiff.jl +https://github.com/JuliaArrays/ArrayInterface.jl +https://github.com/JuliaDiff/FiniteDiff.jl +=# + +function subblockbandrange(A::BandedBlockBandedMatrix) + u, l = subblockbandwidths(A) + return (-l):u +end + +function blockbanded_coloring( + A::Union{BlockBandedMatrix,BandedBlockBandedMatrix}, dim::Integer +) + # consider blocks of columns or rows (let's call them vertices) depending on `dim` + nb_blocks = blocksize(A, dim) + nb_in_block = blocklengths(axes(A, dim)) + first_in_block = blockfirsts(axes(A, dim)) + last_in_block = blocklasts(axes(A, dim)) + color = zeros(Int, size(A, dim)) + + # give a macroscopic color to each block, so that 2 blocks with the same macro color are orthogonal + # same idea as for BandedMatrices + nb_macrocolors = length(blockbandrange(A)) + macrocolor = cycle_range(nb_macrocolors, nb_blocks) + + width = if A isa BandedBlockBandedMatrix + # vertices within a block are colored cleverly using bands + length(subblockbandrange(A)) + else + # vertices within a block are colored naively with distinct micro colors (~ infinite band width) + typemax(Int) + end + + # for each macroscopic color, count how many microscopic colors will be needed + nb_colors_in_macrocolor = zeros(Int, nb_macrocolors) + for mc in 1:nb_macrocolors + largest_nb_in_macrocolor = maximum(nb_in_block[mc:nb_macrocolors:nb_blocks]; init=0) + nb_colors_in_macrocolor[mc] = min(width, largest_nb_in_macrocolor) + end + color_shift_in_macrocolor = vcat(0, cumsum(nb_colors_in_macrocolor)[1:(end - 1)]) + + # assign a microscopic color to each column as a function of its macroscopic color and its position within the block + for b in 1:nb_blocks + block_color = cycle_range(width, nb_in_block[b]) + shift = color_shift_in_macrocolor[macrocolor[b]] + color[first_in_block[b]:last_in_block[b]] .= shift .+ block_color + end + + return color +end + +function SMC.coloring( + A::Union{BlockBandedMatrix,BandedBlockBandedMatrix}, + ::ColoringProblem{:nonsymmetric,:column}, + ::StructuredColoringAlgorithm; + kwargs..., +) + color = blockbanded_coloring(A, 2) + bg = BipartiteGraph(A) + return ColumnColoringResult(A, bg, color) +end + +function SMC.coloring( + A::Union{BlockBandedMatrix,BandedBlockBandedMatrix}, + ::ColoringProblem{:nonsymmetric,:row}, + ::StructuredColoringAlgorithm; + kwargs..., +) + color = blockbanded_coloring(A, 1) + bg = BipartiteGraph(A) + return RowColoringResult(A, bg, color) +end + +end diff --git a/ext/SparseMatrixColoringsCUDAExt.jl b/ext/SparseMatrixColoringsCUDAExt.jl index 2750964d..ed33eece 100644 --- a/ext/SparseMatrixColoringsCUDAExt.jl +++ b/ext/SparseMatrixColoringsCUDAExt.jl @@ -3,23 +3,7 @@ module SparseMatrixColoringsCUDAExt import SparseMatrixColorings as SMC using SparseArrays: SparseMatrixCSC, rowvals, nnz, nzrange using CUDA: CuVector, CuMatrix -using CUDA.CUSPARSE: AbstractCuSparseMatrix, CuSparseMatrixCSC, CuSparseMatrixCSR - -SMC.matrix_versions(A::AbstractCuSparseMatrix) = (A,) - -## Compression (slow, through CPU) - -function SMC.compress( - A::AbstractCuSparseMatrix, result::SMC.AbstractColoringResult{structure,:column} -) where {structure} - return CuMatrix(SMC.compress(SparseMatrixCSC(A), result)) -end - -function SMC.compress( - A::AbstractCuSparseMatrix, result::SMC.AbstractColoringResult{structure,:row} -) where {structure} - return CuMatrix(SMC.compress(SparseMatrixCSC(A), result)) -end +using cuSPARSE: AbstractCuSparseMatrix, CuSparseMatrixCSC, CuSparseMatrixCSR ## CSC Result diff --git a/ext/SparseMatrixColoringsColorsExt.jl b/ext/SparseMatrixColoringsColorsExt.jl index 7f844e89..c934b5e9 100644 --- a/ext/SparseMatrixColoringsColorsExt.jl +++ b/ext/SparseMatrixColoringsColorsExt.jl @@ -271,9 +271,8 @@ function show_colors!( A_ccolor_indices = mod1.(column_colors(res), length(colorscheme)) A_rcolor_indices = mod1.(row_shift .+ row_colors(res), length(colorscheme)) B_ccolor_indices = mod1.(1:maximum(column_colors(res)), length(colorscheme)) - B_rcolor_indices = mod1.( - (row_shift + 1):(row_shift + maximum(row_colors(res))), length(colorscheme) - ) + B_rcolor_indices = + mod1.((row_shift + 1):(row_shift + maximum(row_colors(res))), length(colorscheme)) A_ccolors = colorscheme[A_ccolor_indices] A_rcolors = colorscheme[A_rcolor_indices] B_ccolors = colorscheme[B_ccolor_indices] diff --git a/ext/SparseMatrixColoringsGPUArraysExt.jl b/ext/SparseMatrixColoringsGPUArraysExt.jl new file mode 100644 index 00000000..37574e73 --- /dev/null +++ b/ext/SparseMatrixColoringsGPUArraysExt.jl @@ -0,0 +1,29 @@ +module SparseMatrixColoringsGPUArraysExt + +using GPUArrays: AbstractGPUSparseMatrix, dense_array_type +using SparseArrays: SparseMatrixCSC +import SparseMatrixColorings as SMC + +SMC.matrix_versions(A::AbstractGPUSparseMatrix) = (A,) + +## Compression (slow, through CPU) + +function SMC.compress( + A::AbstractGPUSparseMatrix, result::SMC.AbstractColoringResult{structure,:column} +) where {structure} + A_cpu = SparseMatrixCSC(A) + B_cpu = SMC.compress(A_cpu, result) + B = dense_array_type(A)(B_cpu) + return B +end + +function SMC.compress( + A::AbstractGPUSparseMatrix, result::SMC.AbstractColoringResult{structure,:row} +) where {structure} + A_cpu = SparseMatrixCSC(A) + B_cpu = SMC.compress(A_cpu, result) + B = dense_array_type(A)(B_cpu) + return B +end + +end diff --git a/ext/SparseMatrixColoringsJuMPExt.jl b/ext/SparseMatrixColoringsJuMPExt.jl new file mode 100644 index 00000000..77b74dbd --- /dev/null +++ b/ext/SparseMatrixColoringsJuMPExt.jl @@ -0,0 +1,80 @@ +module SparseMatrixColoringsJuMPExt + +using ADTypes: ADTypes +using JuMP: + Model, + assert_is_solved_and_feasible, + optimize!, + primal_status, + set_silent, + set_start_value, + value, + @variable, + @constraint, + @objective +using JuMP +import MathOptInterface as MOI +using SparseMatrixColorings: + BipartiteGraph, OptimalColoringAlgorithm, nb_vertices, neighbors, pattern, vertices + +function optimal_distance2_coloring( + bg::BipartiteGraph, + ::Val{side}, + optimizer::O; + silent::Bool=true, + assert_solved::Bool=true, +) where {side,O} + other_side = 3 - side + n = nb_vertices(bg, Val(side)) + model = Model(optimizer) + silent && set_silent(model) + # one variable per vertex to color, removing some renumbering symmetries + @variable(model, 1 <= color[i=1:n] <= i, Int) + # one variable to count the number of distinct colors + @variable(model, ncolors, Int) + @constraint(model, [ncolors; color] in MOI.CountDistinct(n + 1)) + # distance-2 coloring: neighbors of the same vertex must have distinct colors + for i in vertices(bg, Val(other_side)) + neigh = neighbors(bg, Val(other_side), i) + @constraint(model, color[neigh] in MOI.AllDifferent(length(neigh))) + end + # minimize the number of distinct colors (can't use maximum because they are not necessarily numbered contiguously) + @objective(model, Min, ncolors) + # actual solving step where time is spent + optimize!(model) + if assert_solved + # assert feasibility and optimality + assert_is_solved_and_feasible(model) + else + # only assert feasibility + @assert primal_status(model) == MOI.FEASIBLE_POINT + end + # native solver solutions are floating point numbers + color_int = round.(Int, value.(color)) + # remap to 1:cmax in case they are not contiguous + true_ncolors = 0 + remap = fill(0, maximum(color_int)) + for c in color_int + if remap[c] == 0 + true_ncolors += 1 + remap[c] = true_ncolors + end + end + return remap[color_int] +end + +function ADTypes.column_coloring(A::AbstractMatrix, algo::OptimalColoringAlgorithm) + bg = BipartiteGraph(A) + return optimal_distance2_coloring( + bg, Val(2), algo.optimizer; algo.silent, algo.assert_solved + ) +end + +function ADTypes.row_coloring(A::AbstractMatrix, algo::OptimalColoringAlgorithm) + bg = BipartiteGraph(A) + return optimal_distance2_coloring( + bg, Val(1), algo.optimizer; algo.silent, algo.assert_solved + ) +end + +end diff --git a/src/SparseMatrixColorings.jl b/src/SparseMatrixColorings.jl index 7c25bac4..74486d80 100644 --- a/src/SparseMatrixColorings.jl +++ b/src/SparseMatrixColorings.jl @@ -14,11 +14,13 @@ using Base.Iterators: Iterators using DocStringExtensions: README, EXPORTS, SIGNATURES, TYPEDEF, TYPEDFIELDS using LinearAlgebra: Adjoint, + Bidiagonal, Diagonal, Hermitian, LowerTriangular, Symmetric, Transpose, + Tridiagonal, UpperTriangular, adjoint, checksquare, @@ -47,15 +49,18 @@ include("graph.jl") include("forest.jl") include("order.jl") include("coloring.jl") +include("postprocessing.jl") include("result.jl") include("matrices.jl") include("interface.jl") include("constant.jl") include("adtypes.jl") include("decompression.jl") +include("structured.jl") include("check.jl") include("examples.jl") include("show_colors.jl") +include("optimal.jl") include("precompile.jl") @@ -63,7 +68,8 @@ export NaturalOrder, RandomOrder, LargestFirst export DynamicDegreeBasedOrder, SmallestLast, IncidenceDegree, DynamicLargestFirst export PerfectEliminationOrder export ColoringProblem, GreedyColoringAlgorithm, AbstractColoringResult -export ConstantColoringAlgorithm +export ConstantColoringAlgorithm, StructuredColoringAlgorithm +export OptimalColoringAlgorithm export coloring, fast_coloring export column_colors, row_colors, ncolors export column_groups, row_groups diff --git a/src/adtypes.jl b/src/adtypes.jl index 7500b1e2..7c464433 100644 --- a/src/adtypes.jl +++ b/src/adtypes.jl @@ -1,21 +1,33 @@ -function coloring( - A::AbstractMatrix, - ::ColoringProblem{:nonsymmetric,:column}, - algo::ADTypes.NoColoringAlgorithm; - kwargs..., -) - bg = BipartiteGraph(A) - color = convert(Vector{eltype(bg)}, ADTypes.column_coloring(A, algo)) - return ColumnColoringResult(A, bg, color) -end +## From ADTypes to SMC function coloring( A::AbstractMatrix, - ::ColoringProblem{:nonsymmetric,:row}, - algo::ADTypes.NoColoringAlgorithm; - kwargs..., -) - bg = BipartiteGraph(A) - color = convert(Vector{eltype(bg)}, ADTypes.row_coloring(A, algo)) - return RowColoringResult(A, bg, color) + problem::ColoringProblem{structure,partition}, + algo::ADTypes.AbstractColoringAlgorithm; + decompression_eltype::Type{R}=Float64, + symmetric_pattern::Bool=false, +) where {structure,partition,R} + symmetric_pattern = symmetric_pattern || A isa Union{Symmetric,Hermitian} + if structure == :nonsymmetric + if partition == :column + forced_colors = ADTypes.column_coloring(A, algo) + elseif partition == :row + forced_colors = ADTypes.row_coloring(A, algo) + else + # TODO: improve once https://github.com/SciML/ADTypes.jl/issues/69 is done + A_and_Aᵀ, _ = bidirectional_pattern(A; symmetric_pattern) + forced_colors = ADTypes.symmetric_coloring(A_and_Aᵀ, algo) + end + else + forced_colors = ADTypes.symmetric_coloring(A, algo) + end + return _coloring( + WithResult(), + A, + problem, + GreedyColoringAlgorithm(), + R, + symmetric_pattern; + forced_colors, + ) end diff --git a/src/check.jl b/src/check.jl index f47f0958..db5c0a34 100644 --- a/src/check.jl +++ b/src/check.jl @@ -86,11 +86,15 @@ A partition of the columns of a symmetric matrix `A` is _symmetrically orthogona 1. the group containing the column `A[:, j]` has no other column with a nonzero in row `i` 2. the group containing the column `A[:, i]` has no other column with a nonzero in row `j` +It is equivalent to a __star coloring__. + !!! warning This function is not coded with efficiency in mind, it is designed for small-scale tests. # References +> [_On the Estimation of Sparse Hessian Matrices_](https://doi.org/10.1137/0716078), Powell and Toint (1979) +> [_Estimation of sparse hessian matrices and graph coloring problems_](https://doi.org/10.1007/BF02612334), Coleman and Moré (1984) > [_What Color Is Your Jacobian? Graph Coloring for Computing Derivatives_](https://epubs.siam.org/doi/10.1137/S0036144504444711), Gebremedhin et al. (2005) """ function symmetrically_orthogonal_columns( @@ -102,7 +106,7 @@ function symmetrically_orthogonal_columns( end issymmetric(A) || return false group = group_by_color(color) - for i in axes(A, 2), j in axes(A, 2) + for i in axes(A, 1), j in axes(A, 2) iszero(A[i, j]) && continue ci, cj = color[i], color[j] check = _bilateral_check( @@ -126,6 +130,8 @@ A bipartition of the rows and columns of a matrix `A` is _structurally biorthogo 1. the group containing the column `A[:, j]` has no other column with a nonzero in row `i` 2. the group containing the row `A[i, :]` has no other row with a nonzero in column `j` +It is equivalent to a __star bicoloring__. + !!! warning This function is not coded with efficiency in mind, it is designed for small-scale tests. """ @@ -138,8 +144,8 @@ function structurally_biorthogonal( if !proper_length_bicoloring(A, row_color, column_color; verbose) return false end - column_group = group_by_color(column_color) row_group = group_by_color(row_color) + column_group = group_by_color(column_color) for i in axes(A, 1), j in axes(A, 2) iszero(A[i, j]) && continue ci, cj = row_color[i], column_color[j] @@ -261,6 +267,174 @@ function directly_recoverable_columns( return true end +""" + substitutable_columns( + A::AbstractMatrix, rank_nonzeros::AbstractMatrix, color::AbstractVector{<:Integer}; + verbose=false + ) + +Return `true` if coloring the columns of the symmetric matrix `A` with the vector `color` results in a partition that is substitutable, and `false` otherwise. +For all nonzeros `A[i, j]`, `rank_nonzeros[i, j]` provides its rank of recovery. + +A partition of the columns of a symmetric matrix `A` is _substitutable_ if, for every nonzero element `A[i, j]`, either of the following statements holds: + +1. the group containing the column `A[:, j]` has all nonzeros in row `i` ordered before `A[i, j]` +2. the group containing the column `A[:, i]` has all nonzeros in row `j` ordered before `A[i, j]` + +It is equivalent to an __acyclic coloring__. + +!!! warning + This function is not coded with efficiency in mind, it is designed for small-scale tests. + +# References + +> [_On the Estimation of Sparse Hessian Matrices_](https://doi.org/10.1137/0716078), Powell and Toint (1979) +> [_The Cyclic Coloring Problem and Estimation of Sparse Hessian Matrices_](https://doi.org/10.1137/0607026), Coleman and Cai (1986) +> [_What Color Is Your Jacobian? Graph Coloring for Computing Derivatives_](https://epubs.siam.org/doi/10.1137/S0036144504444711), Gebremedhin et al. (2005) +""" +function substitutable_columns( + A::AbstractMatrix, + rank_nonzeros::AbstractMatrix, + color::AbstractVector{<:Integer}; + verbose::Bool=false, +) + checksquare(A) + if !proper_length_coloring(A, color; verbose) + return false + end + issymmetric(A) || return false + group = group_by_color(color) + for i in axes(A, 1), j in axes(A, 2) + iszero(A[i, j]) && continue + ci, cj = color[i], color[j] + check = _substitutable_check( + A, rank_nonzeros; i, j, ci, cj, row_group=group, column_group=group, verbose + ) + !check && return false + end + return true +end + +""" + substitutable_bidirectional( + A::AbstractMatrix, rank_nonzeros::AbstractMatrix, row_color::AbstractVector{<:Integer}, column_color::AbstractVector{<:Integer}; + verbose=false + ) + +Return `true` if bicoloring of the matrix `A` with the vectors `row_color` and `column_color` results in a bipartition that is substitutable, and `false` otherwise. +For all nonzeros `A[i, j]`, `rank_nonzeros[i, j]` provides its rank of recovery. + +A bipartition of the rows and columns of a matrix `A` is _substitutable_ if, for every nonzero element `A[i, j]`, either of the following statements holds: + +1. the group containing the column `A[:, j]` has all nonzeros in row `i` ordered before `A[i, j]` +2. the group containing the row `A[i, :]` has all nonzeros in column `j` ordered before `A[i, j]` + +It is equivalent to an __acyclic bicoloring__. + +!!! warning + This function is not coded with efficiency in mind, it is designed for small-scale tests. +""" +function substitutable_bidirectional( + A::AbstractMatrix, + rank_nonzeros::AbstractMatrix, + row_color::AbstractVector{<:Integer}, + column_color::AbstractVector{<:Integer}; + verbose::Bool=false, +) + if !proper_length_bicoloring(A, row_color, column_color; verbose) + return false + end + row_group = group_by_color(row_color) + column_group = group_by_color(column_color) + for i in axes(A, 1), j in axes(A, 2) + iszero(A[i, j]) && continue + ci, cj = row_color[i], column_color[j] + check = _substitutable_check( + A, rank_nonzeros; i, j, ci, cj, row_group, column_group, verbose + ) + !check && return false + end + return true +end + +function _substitutable_check( + A::AbstractMatrix, + rank_nonzeros::AbstractMatrix; + i::Integer, + j::Integer, + ci::Integer, + cj::Integer, + row_group::AbstractVector, + column_group::AbstractVector, + verbose::Bool, +) + order_ij = rank_nonzeros[i, j] + k_row = 0 + k_column = 0 + if ci != 0 + for k in row_group[ci] + (k == i) && continue + if !iszero(A[k, j]) + order_kj = rank_nonzeros[k, j] + @assert !iszero(order_kj) + if order_kj > order_ij + k_row = k + end + end + end + end + if cj != 0 + for k in column_group[cj] + (k == j) && continue + if !iszero(A[i, k]) + order_ik = rank_nonzeros[i, k] + @assert !iszero(order_ik) + if order_ik > order_ij + k_column = k + end + end + end + end + if ci == 0 && cj == 0 + if verbose + @warn """ + For coefficient (i=$i, j=$j) with colors (ci=$ci, cj=$cj): + - Row color ci=$ci is neutral. + - Column color cj=$cj is neutral. + """ + end + return false + elseif ci == 0 && !iszero(k_column) + if verbose + @warn """ + For coefficient (i=$i, j=$j) with colors (ci=$ci, cj=$cj): + - Row color ci=$ci is neutral. + - For the column $k_column in column color cj=$cj, A[$i, $k_column] is ordered after A[$i, $j]. + """ + end + return false + elseif cj == 0 && !iszero(k_row) + if verbose + @warn """ + For coefficient (i=$i, j=$j) with colors (ci=$ci, cj=$cj): + - For the row $k_row in row color ci=$ci, A[$k_row, $j] is ordered after A[$i, $j]. + - Column color cj=$cj is neutral. + """ + end + return false + elseif !iszero(k_row) && !iszero(k_column) + if verbose + @warn """ + For coefficient (i=$i, j=$j) with colors (ci=$ci, cj=$cj): + - For the row $k_row in row color ci=$ci, A[$k_row, $j] is ordered after A[$i, $j]. + - For the column $k_column in column color cj=$cj, A[$i, $k_column] is ordered after A[$i, $j]. + """ + end + return false + end + return true +end + """ valid_dynamic_order(g::AdjacencyGraph, π::AbstractVector{<:Integer}, order::DynamicDegreeBasedOrder) valid_dynamic_order(bg::AdjacencyGraph, ::Val{side}, π::AbstractVector{<:Integer}, order::DynamicDegreeBasedOrder) @@ -326,3 +500,75 @@ function valid_dynamic_order( end return true end + +""" + rank_nonzeros_from_trees(result::TreeSetColoringResult) + rank_nonzeros_from_trees(result::BicoloringResult) + +Construct a sparse matrix `rank_nonzeros` that assigns a unique recovery rank +to each nonzero coefficient associated with an acyclic coloring or bicoloring. + +For every nonzero entry `result.A[i, j]`, `rank_nonzeros[i, j]` stores a strictly positive +integer representing the order in which this coefficient is recovered during the decompression. +A larger value means the coefficient is recovered later. + +This ranking is used to test substitutability (acyclicity) of colorings: +for a given nonzero `result.A[i, j]`, the ranks allow one to check whether all competing +nonzeros in the same row or column (within a color group) are recovered before it. +""" +function rank_nonzeros_from_trees end + +function rank_nonzeros_from_trees(result::TreeSetColoringResult) + (; A, ag, reverse_bfs_orders, diagonal_indices, tree_edge_indices, nt) = result + (; S) = ag + m, n = size(A) + nnzS = nnz(S) + nzval = zeros(Int, nnzS) + rank_nonzeros = SparseMatrixCSC(n, n, S.colptr, S.rowval, nzval) + counter = 0 + for i in diagonal_indices + counter += 1 + rank_nonzeros[i, i] = counter + end + for k in 1:nt + first = tree_edge_indices[k] + last = tree_edge_indices[k + 1] - 1 + for pos in first:last + (i, j) = reverse_bfs_orders[pos] + counter += 1 + rank_nonzeros[i, j] = counter + rank_nonzeros[j, i] = counter + end + end + return rank_nonzeros +end + +function rank_nonzeros_from_trees(result::BicoloringResult) + (; A, abg, row_color, column_color, symmetric_result, large_colptr, large_rowval) = + result + @assert symmetric_result isa TreeSetColoringResult + (; ag, reverse_bfs_orders, tree_edge_indices, nt) = symmetric_result + (; S) = ag + m, n = size(A) + nnzA = nnz(S) ÷ 2 + nzval = zeros(Int, nnzA) + colptr = large_colptr[1:(n + 1)] + rowval = large_rowval[1:nnzA] + rowval .-= n + rank_nonzeros = SparseMatrixCSC(m, n, colptr, rowval, nzval) + counter = 0 + for k in 1:nt + first = tree_edge_indices[k] + last = tree_edge_indices[k + 1] - 1 + for pos in first:last + (i, j) = reverse_bfs_orders[pos] + counter += 1 + if i > j + rank_nonzeros[i - n, j] = counter + else + rank_nonzeros[j - n, i] = counter + end + end + end + return rank_nonzeros +end diff --git a/src/coloring.jl b/src/coloring.jl index 7eec55c0..699dbc64 100644 --- a/src/coloring.jl +++ b/src/coloring.jl @@ -1,5 +1,10 @@ +struct InvalidColoringError <: Exception end + """ - partial_distance2_coloring(bg::BipartiteGraph, ::Val{side}, vertices_in_order::AbstractVector) + partial_distance2_coloring( + bg::BipartiteGraph, ::Val{side}, vertices_in_order::AbstractVector; + forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing + ) Compute a distance-2 coloring of the given `side` (`1` or `2`) in the bipartite graph `bg` and return a vector of integer colors. @@ -7,6 +12,8 @@ A _distance-2 coloring_ is such that two vertices have different colors if they The vertices are colored in a greedy fashion, following the order supplied. +The optional `forced_colors` keyword argument is used to enforce predefined vertex colors (e.g. coming from another optimization algorithm) but still run the distance-2 coloring procedure to verify correctness. + # See also - [`BipartiteGraph`](@ref) @@ -17,11 +24,16 @@ The vertices are colored in a greedy fashion, following the order supplied. > [_What Color Is Your Jacobian? Graph Coloring for Computing Derivatives_](https://epubs.siam.org/doi/10.1137/S0036144504444711), Gebremedhin et al. (2005), Algorithm 3.2 """ function partial_distance2_coloring( - bg::BipartiteGraph{T}, ::Val{side}, vertices_in_order::AbstractVector{<:Integer} + bg::BipartiteGraph{T}, + ::Val{side}, + vertices_in_order::AbstractVector{<:Integer}; + forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing, ) where {T,side} color = Vector{T}(undef, nb_vertices(bg, Val(side))) forbidden_colors = Vector{T}(undef, nb_vertices(bg, Val(side))) - partial_distance2_coloring!(color, forbidden_colors, bg, Val(side), vertices_in_order) + partial_distance2_coloring!( + color, forbidden_colors, bg, Val(side), vertices_in_order; forced_colors + ) return color end @@ -30,7 +42,8 @@ function partial_distance2_coloring!( forbidden_colors::AbstractVector{<:Integer}, bg::BipartiteGraph, ::Val{side}, - vertices_in_order::AbstractVector{<:Integer}, + vertices_in_order::AbstractVector{<:Integer}; + forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing, ) where {side} color .= 0 forbidden_colors .= 0 @@ -44,17 +57,32 @@ function partial_distance2_coloring!( end end end - for i in eachindex(forbidden_colors) - if forbidden_colors[i] != v - color[v] = i - break + if isnothing(forced_colors) + for i in eachindex(forbidden_colors) + if forbidden_colors[i] != v + color[v] = i + break + end + end + else + f = forced_colors[v] + if ( + (f == 0 && length(neighbors(bg, Val(side), v)) > 0) || + (f > 0 && forbidden_colors[f] == v) + ) + throw(InvalidColoringError()) + else + color[v] = f end end end end """ - star_coloring(g::AdjacencyGraph, vertices_in_order::AbstractVector, postprocessing::Bool) + star_coloring( + g::AdjacencyGraph, vertices_in_order::AbstractVector, postprocessing::Bool; + forced_colors::Union{AbstractVector,Nothing}=nothing + ) Compute a star coloring of all vertices in the adjacency graph `g` and return a tuple `(color, star_set)`, where @@ -67,6 +95,8 @@ The vertices are colored in a greedy fashion, following the order supplied. If `postprocessing=true`, some colors might be replaced with `0` (the "neutral" color) as long as they are not needed during decompression. +The optional `forced_colors` keyword argument is used to enforce predefined vertex colors (e.g. coming from another optimization algorithm) but still run the star coloring procedure to verify correctness and build auxiliary data structures, useful during decompression. + # See also - [`AdjacencyGraph`](@ref) @@ -77,7 +107,10 @@ If `postprocessing=true`, some colors might be replaced with `0` (the "neutral" > [_New Acyclic and Star Coloring Algorithms with Application to Computing Hessians_](https://epubs.siam.org/doi/abs/10.1137/050639879), Gebremedhin et al. (2007), Algorithm 4.1 """ function star_coloring( - g::AdjacencyGraph{T}, vertices_in_order::AbstractVector{<:Integer}, postprocessing::Bool + g::AdjacencyGraph{T}, + vertices_in_order::AbstractVector{<:Integer}, + postprocessing::Bool; + forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing, ) where {T<:Integer} # Initialize data structures nv = nb_vertices(g) @@ -91,7 +124,7 @@ function star_coloring( for v in vertices_in_order for (w, index_vw) in neighbors_with_edge_indices(g, v) - !has_diagonal(g) || (v == w && continue) + augmented_graph(g) || (v == w && continue) color_w = color[w] iszero(color_w) && continue forbidden_colors[color_w] = v @@ -106,7 +139,7 @@ function star_coloring( else first_neighbor[color[w]] = (v, w, index_vw) for (x, index_wx) in neighbors_with_edge_indices(g, w) - !has_diagonal(g) || (w == x && continue) + augmented_graph(g) || (w == x && continue) color_x = color[x] (x == v || iszero(color_x)) && continue if x == hub[star[index_wx]] # potential Case 2 (which is always false for trivial stars with two vertices, since the associated hub is negative) @@ -115,10 +148,18 @@ function star_coloring( end end end - for i in eachindex(forbidden_colors) - if forbidden_colors[i] != v - color[v] = i - break + if isnothing(forced_colors) + for i in eachindex(forbidden_colors) + if forbidden_colors[i] != v + color[v] = i + break + end + end + else + if forbidden_colors[forced_colors[v]] == v # TODO: handle forced_colors[v] == 0 + throw(InvalidColoringError()) + else + color[v] = forced_colors[v] end end _update_stars!(star, hub, g, v, color, first_neighbor) @@ -143,7 +184,7 @@ function _treat!( color::AbstractVector{<:Integer}, ) for x in neighbors(g, w) - !has_diagonal(g) || (w == x && continue) + augmented_graph(g) || (w == x && continue) color_x = color[x] iszero(color_x) && continue forbidden_colors[color_x] = v @@ -163,12 +204,12 @@ function _update_stars!( first_neighbor::AbstractVector{<:Tuple}, ) for (w, index_vw) in neighbors_with_edge_indices(g, v) - !has_diagonal(g) || (v == w && continue) + augmented_graph(g) || (v == w && continue) color_w = color[w] iszero(color_w) && continue x_exists = false for (x, index_wx) in neighbors_with_edge_indices(g, w) - !has_diagonal(g) || (w == x && continue) + augmented_graph(g) || (w == x && continue) if x != v && color[x] == color[v] # vw, wx ∈ E star_wx = star[index_wx] hub[star_wx] = w # this may already be true @@ -245,23 +286,22 @@ function acyclic_coloring( for v in vertices_in_order for w in neighbors(g, v) - !has_diagonal(g) || (v == w && continue) + augmented_graph(g) || (v == w && continue) color_w = color[w] iszero(color_w) && continue forbidden_colors[color_w] = v end for w in neighbors(g, v) - !has_diagonal(g) || (v == w && continue) + augmented_graph(g) || (v == w && continue) iszero(color[w]) && continue for (x, index_wx) in neighbors_with_edge_indices(g, w) - !has_diagonal(g) || (w == x && continue) + augmented_graph(g) || (w == x && continue) color_x = color[x] iszero(color_x) && continue if forbidden_colors[color_x] != v _prevent_cycle!( v, w, - x, index_wx, color_x, first_visit_to_tree, @@ -271,6 +311,7 @@ function acyclic_coloring( end end end + # TODO: handle forced colors for i in eachindex(forbidden_colors) if forbidden_colors[i] != v color[v] = i @@ -278,20 +319,20 @@ function acyclic_coloring( end end for (w, index_vw) in neighbors_with_edge_indices(g, v) # grow two-colored stars around the vertex v - !has_diagonal(g) || (v == w && continue) + augmented_graph(g) || (v == w && continue) color_w = color[w] iszero(color_w) && continue _grow_star!(v, w, index_vw, color_w, first_neighbor, forest) end for (w, index_vw) in neighbors_with_edge_indices(g, v) - !has_diagonal(g) || (v == w && continue) + augmented_graph(g) || (v == w && continue) iszero(color[w]) && continue for (x, index_wx) in neighbors_with_edge_indices(g, w) - !has_diagonal(g) || (w == x && continue) + augmented_graph(g) || (w == x && continue) color_x = color[x] (x == v || iszero(color_x)) && continue if color_x == color[v] - _merge_trees!(v, w, x, index_vw, index_wx, forest) # merge trees T₁ ∋ vw and T₂ ∋ wx if T₁ != T₂ + _merge_trees!(index_vw, index_wx, forest) # merge trees T₁ ∋ vw and T₂ ∋ wx if T₁ != T₂ end end end @@ -312,7 +353,6 @@ function _prevent_cycle!( # not modified v::Integer, w::Integer, - x::Integer, index_wx::Integer, color_x::Integer, # modified @@ -354,9 +394,6 @@ end function _merge_trees!( # not modified - v::Integer, - w::Integer, - x::Integer, index_vw::Integer, index_wx::Integer, # modified @@ -380,9 +417,23 @@ Encode a set of 2-colored trees resulting from the [`acyclic_coloring`](@ref) al $TYPEDFIELDS """ struct TreeSet{T} + """ + For a tree index `1 <= k <= nt`, the list + `list = reverse_bfs_order[tree_edge_indices[k]:(tree_edge_indices[k+1]-1)]` is a list of edges + `list[i] = (leaf, inner)` of the `k`th tree such that `leaf` is a leaf of the tree containing + the edges `list[i:end]`. + From an other point of view, `reverse(list)` contains the edges in the order of a breadth first + (BFS) traversal order of the `k`th tree starting from a depth-minimizing root. + """ reverse_bfs_orders::Vector{Tuple{T,T}} + "For a tree index `1 <= k <= nt`, `is_star[k]` indicates whether the `k`th three is a star." is_star::Vector{Bool} + """ + `tree_edge_indices[1]` is one and `tree_edge_indices[k+1] - tree_edge_indices[k]` is the number of edges in the `k`th tree. + One can think of it as a kind of fused vector of offsets (similar to the `colptr` field of `SparseMatrixCSC`) of all trees together. + """ tree_edge_indices::Vector{T} + "numbers of 2-colored trees for which trees sharing the same 2 colors have disjoint edges" nt::T end @@ -390,6 +441,7 @@ function TreeSet( g::AdjacencyGraph{T}, forest::Forest{T}, buffer::AbstractVector{T}, + # The value of `reverse_bfs_orders` is ignored, we just provide the storage for it (or reuse memory allocated during acyclic coloring) reverse_bfs_orders::Vector{Tuple{T,T}}, ne::Integer, ) where {T} @@ -524,7 +576,15 @@ function TreeSet( # Number of edges treated num_edges_treated = zero(T) - # reverse_bfs_orders contains the reverse breadth first (BFS) traversal order for each tree in the forest + # The `rank` of the `k`th tree encoded in `forest` does not correspond + # to the depth of the tree rooted as the root encoded in `forest` because + # `forest.parents[u] = v` only needs a path to exists from `u` to `v` but + # there may not be an edge `(u, v)`. + # We also want a root `r` that minimizes the depth of the tree rooted at + # `r`. To achieve this, we start from each leaf and remove the corresponding + # edges. We then look at all leaves of the corresponding graphs and repeat + # the process until there is only one vertex left. This vertex will then be + # a depth-minimizing root. for k in 1:nt # Initialize the queue to store the leaves queue_start = 1 @@ -601,150 +661,3 @@ function TreeSet( return TreeSet(reverse_bfs_orders, is_star, tree_edge_indices, nt) end - -## Postprocessing, mirrors decompression code - -function postprocess!( - color::AbstractVector{<:Integer}, - star_or_tree_set::Union{StarSet,TreeSet}, - g::AdjacencyGraph, - offsets::AbstractVector{<:Integer}, -) - S = pattern(g) - edge_to_index = edge_indices(g) - # flag which colors are actually used during decompression - nb_colors = maximum(color) - color_used = zeros(Bool, nb_colors) - - # nonzero diagonal coefficients force the use of their respective color (there can be no neutral colors if the diagonal is fully nonzero) - if has_diagonal(g) - for i in axes(S, 1) - if !iszero(S[i, i]) - color_used[color[i]] = true - end - end - end - - if star_or_tree_set isa StarSet - # only the colors of the hubs are used - (; star, hub) = star_or_tree_set - nb_trivial_stars = 0 - - # Iterate through all non-trivial stars - for s in eachindex(hub) - h = hub[s] - if h > 0 - color_used[color[h]] = true - else - nb_trivial_stars += 1 - end - end - - # Process the trivial stars (if any) - if nb_trivial_stars > 0 - rvS = rowvals(S) - for j in axes(S, 2) - for k in nzrange(S, j) - i = rvS[k] - if i > j - index_ij = edge_to_index[k] - s = star[index_ij] - h = hub[s] - if h < 0 - h = abs(h) - spoke = h == j ? i : j - if color_used[color[spoke]] - # Switch the hub and the spoke to possibly avoid adding one more used color - hub[s] = spoke - else - # Keep the current hub - color_used[color[h]] = true - end - end - end - end - end - end - else - # only the colors of non-leaf vertices are used - (; reverse_bfs_orders, is_star, tree_edge_indices, nt) = star_or_tree_set - nb_trivial_trees = 0 - - # Iterate through all non-trivial trees - for k in 1:nt - # Position of the first edge in the tree - first = tree_edge_indices[k] - - # Total number of edges in the tree - ne_tree = tree_edge_indices[k + 1] - first - - # Check if we have more than one edge in the tree (non-trivial tree) - if ne_tree > 1 - # Determine if the tree is a star - if is_star[k] - # It is a non-trivial star and only the color of the hub is needed - (_, hub) = reverse_bfs_orders[first] - color_used[color[hub]] = true - else - # It is not a star and both colors are needed during the decompression - (i, j) = reverse_bfs_orders[first] - color_used[color[i]] = true - color_used[color[j]] = true - end - else - nb_trivial_trees += 1 - end - end - - # Process the trivial trees (if any) - if nb_trivial_trees > 0 - for k in 1:nt - # Position of the first edge in the tree - first = tree_edge_indices[k] - - # Total number of edges in the tree - ne_tree = tree_edge_indices[k + 1] - first - - # Check if we have exactly one edge in the tree - if ne_tree == 1 - (i, j) = reverse_bfs_orders[first] - if color_used[color[i]] - # Make i the root to avoid possibly adding one more used color - # Switch it with the (only) leaf - reverse_bfs_orders[first] = (j, i) - else - # Keep j as the root - color_used[color[j]] = true - end - end - end - end - end - - # if at least one of the colors is useless, modify the color assignments of vertices - if any(!, color_used) - num_colors_useless = 0 - - # determine what are the useless colors and compute the offsets - for ci in 1:nb_colors - if color_used[ci] - offsets[ci] = num_colors_useless - else - num_colors_useless += 1 - end - end - - # assign the neutral color to every vertex with a useless color and remap the colors - for i in eachindex(color) - ci = color[i] - if !color_used[ci] - # assign the neutral color - color[i] = 0 - else - # remap the color to not have any gap - color[i] -= offsets[ci] - end - end - end - return color -end diff --git a/src/constant.jl b/src/constant.jl index d79caab0..56bf4101 100644 --- a/src/constant.jl +++ b/src/constant.jl @@ -4,20 +4,24 @@ Coloring algorithm which always returns the same precomputed vector of colors. Useful when the optimal coloring of a matrix can be determined a priori due to its specific structure (e.g. banded). -It is passed as an argument to the main function [`coloring`](@ref), but will only work if the associated `problem` has `:nonsymmetric` structure. -Indeed, for symmetric coloring problems, we need more than just the vector of colors to allow fast decompression. +It is passed as an argument to the main function [`coloring`](@ref), but will only work if the associated `problem` has a `:column` or `:row` partition. # Constructors ConstantColoringAlgorithm{partition}(matrix_template, color) - ConstantColoringAlgorithm(matrix_template, color; partition=:column) + ConstantColoringAlgorithm{partition,structure}(matrix_template, color) + ConstantColoringAlgorithm( + matrix_template, color; + structure=:nonsymmetric, partition=:column + ) - `partition::Symbol`: either `:row` or `:column`. +- `structure::Symbol`: either `:nonsymmetric` or `:symmetric`. - `matrix_template::AbstractMatrix`: matrix for which the vector of colors was precomputed (the algorithm will only accept matrices of the exact same size). - `color::Vector{<:Integer}`: vector of integer colors, one for each row or column (depending on `partition`). !!! warning - The second constructor (based on keyword arguments) is type-unstable. + The constructor based on keyword arguments is type-unstable if these arguments are not compile-time constants. We do not necessarily verify consistency between the matrix template and the vector of colors, this is the responsibility of the user. @@ -63,71 +67,68 @@ julia> column_colors(result) - [`ADTypes.column_coloring`](@extref ADTypes.column_coloring) - [`ADTypes.row_coloring`](@extref ADTypes.row_coloring) +- [`ADTypes.symmetric_coloring`](@extref ADTypes.symmetric_coloring) """ -struct ConstantColoringAlgorithm{ - partition, - M<:AbstractMatrix, - T<:Integer, - R<:AbstractColoringResult{:nonsymmetric,partition,:direct}, -} <: ADTypes.AbstractColoringAlgorithm +struct ConstantColoringAlgorithm{partition,structure,M<:AbstractMatrix,T<:Integer} <: + ADTypes.AbstractColoringAlgorithm matrix_template::M color::Vector{T} - result::R -end -function ConstantColoringAlgorithm{:column}( - matrix_template::AbstractMatrix, color::Vector{<:Integer} -) - bg = BipartiteGraph(matrix_template) - result = ColumnColoringResult(matrix_template, bg, color) - T, M, R = eltype(bg), typeof(matrix_template), typeof(result) - return ConstantColoringAlgorithm{:column,M,T,R}(matrix_template, color, result) + function ConstantColoringAlgorithm{partition,structure}( + matrix_template::AbstractMatrix, color::Vector{<:Integer} + ) where {partition,structure} + check_valid_problem(structure, partition) + return new{partition,structure,typeof(matrix_template),eltype(color)}( + matrix_template, color + ) + end end -function ConstantColoringAlgorithm{:row}( +function ConstantColoringAlgorithm{partition}( matrix_template::AbstractMatrix, color::Vector{<:Integer} -) - bg = BipartiteGraph(matrix_template) - result = RowColoringResult(matrix_template, bg, color) - T, M, R = eltype(bg), typeof(matrix_template), typeof(result) - return ConstantColoringAlgorithm{:row,M,T,R}(matrix_template, color, result) +) where {partition} + return ConstantColoringAlgorithm{partition,:nonsymmetric}(matrix_template, color) end function ConstantColoringAlgorithm( - matrix_template::AbstractMatrix, color::Vector{<:Integer}; partition::Symbol=:column + matrix_template::AbstractMatrix, + color::Vector{<:Integer}; + structure::Symbol=:nonsymmetric, + partition::Symbol=:column, ) - return ConstantColoringAlgorithm{partition}(matrix_template, color) + return ConstantColoringAlgorithm{partition,structure}(matrix_template, color) end -function coloring( - A::AbstractMatrix, - ::ColoringProblem{:nonsymmetric,partition}, - algo::ConstantColoringAlgorithm{partition}; - decompression_eltype::Type=Float64, - symmetric_pattern::Bool=false, -) where {partition} - (; matrix_template, result) = algo +function check_template(algo::ConstantColoringAlgorithm, A::AbstractMatrix) + (; matrix_template) = algo if size(A) != size(matrix_template) throw( DimensionMismatch( "`ConstantColoringAlgorithm` expected matrix of size $(size(matrix_template)) but got matrix of size $(size(A))", ), ) - else - return result end end function ADTypes.column_coloring( - A::AbstractMatrix, algo::ConstantColoringAlgorithm{:column} + A::AbstractMatrix, algo::ConstantColoringAlgorithm{:column,:nonsymmetric} +) + check_template(algo, A) + return algo.color +end + +function ADTypes.row_coloring( + A::AbstractMatrix, algo::ConstantColoringAlgorithm{:row,:nonsymmetric} ) - problem = ColoringProblem{:nonsymmetric,:column}() - result = coloring(A, problem, algo) - return column_colors(result) + check_template(algo, A) + return algo.color end -function ADTypes.row_coloring(A::AbstractMatrix, algo::ConstantColoringAlgorithm) - problem = ColoringProblem{:nonsymmetric,:row}() - result = coloring(A, problem, algo) - return row_colors(result) +function ADTypes.symmetric_coloring( + A::AbstractMatrix, algo::ConstantColoringAlgorithm{:column,:symmetric} +) + check_template(algo, A) + return algo.color end + +# TODO: handle bidirectional once https://github.com/SciML/ADTypes.jl/issues/69 is done diff --git a/src/decompression.jl b/src/decompression.jl index a3dd471a..849f0f6b 100644 --- a/src/decompression.jl +++ b/src/decompression.jl @@ -106,7 +106,7 @@ function compress( return Br, Bc end -struct UnsupportedDecompressionError +struct UnsupportedDecompressionError <: Exception msg::String end @@ -339,9 +339,9 @@ end ## ColumnColoringResult function decompress!(A::AbstractMatrix, B::AbstractMatrix, result::ColumnColoringResult) - (; color) = result - S = result.bg.S2 - check_same_pattern(A, S) + (; bg, color) = result + S = bg.S2 + check_compatible_pattern(A, bg) fill!(A, zero(eltype(A))) rvS = rowvals(S) for j in axes(S, 2) @@ -357,9 +357,9 @@ end function decompress_single_color!( A::AbstractMatrix, b::AbstractVector, c::Integer, result::ColumnColoringResult ) - (; group) = result - S = result.bg.S2 - check_same_pattern(A, S) + (; bg, group) = result + S = bg.S2 + check_compatible_pattern(A, bg) rvS = rowvals(S) for j in group[c] for k in nzrange(S, j) @@ -371,9 +371,9 @@ function decompress_single_color!( end function decompress!(A::SparseMatrixCSC, B::AbstractMatrix, result::ColumnColoringResult) - (; compressed_indices) = result - S = result.bg.S2 - check_same_pattern(A, S) + (; bg, compressed_indices) = result + S = bg.S2 + check_compatible_pattern(A, bg) nzA = nonzeros(A) for k in eachindex(nzA, compressed_indices) nzA[k] = B[compressed_indices[k]] @@ -384,9 +384,9 @@ end function decompress_single_color!( A::SparseMatrixCSC, b::AbstractVector, c::Integer, result::ColumnColoringResult ) - (; group) = result - S = result.bg.S2 - check_same_pattern(A, S) + (; bg, group) = result + S = bg.S2 + check_compatible_pattern(A, bg) rvS = rowvals(S) nzA = nonzeros(A) for j in group[c] @@ -401,9 +401,9 @@ end ## RowColoringResult function decompress!(A::AbstractMatrix, B::AbstractMatrix, result::RowColoringResult) - (; color) = result - S = result.bg.S2 - check_same_pattern(A, S) + (; bg, color) = result + S = bg.S2 + check_compatible_pattern(A, bg) fill!(A, zero(eltype(A))) rvS = rowvals(S) for j in axes(S, 2) @@ -419,9 +419,9 @@ end function decompress_single_color!( A::AbstractMatrix, b::AbstractVector, c::Integer, result::RowColoringResult ) - (; group) = result - S, Sᵀ = result.bg.S2, result.bg.S1 - check_same_pattern(A, S) + (; bg, group) = result + S, Sᵀ = bg.S2, bg.S1 + check_compatible_pattern(A, bg) rvSᵀ = rowvals(Sᵀ) for i in group[c] for k in nzrange(Sᵀ, i) @@ -433,9 +433,9 @@ function decompress_single_color!( end function decompress!(A::SparseMatrixCSC, B::AbstractMatrix, result::RowColoringResult) - (; compressed_indices) = result - S = result.bg.S2 - check_same_pattern(A, S) + (; bg, compressed_indices) = result + S = bg.S2 + check_compatible_pattern(A, bg) nzA = nonzeros(A) for k in eachindex(nzA, compressed_indices) nzA[k] = B[compressed_indices[k]] @@ -450,7 +450,7 @@ function decompress!( ) (; ag, compressed_indices) = result (; S) = ag - uplo == :F && check_same_pattern(A, S) + check_compatible_pattern(A, ag, uplo) fill!(A, zero(eltype(A))) rvS = rowvals(S) @@ -474,7 +474,7 @@ function decompress_single_color!( ) (; ag, compressed_indices, group) = result (; S) = ag - uplo == :F && check_same_pattern(A, S) + check_compatible_pattern(A, ag, uplo) lower_index = (c - 1) * S.n + 1 upper_index = c * S.n @@ -508,8 +508,8 @@ function decompress!( (; ag, compressed_indices) = result (; S) = ag nzA = nonzeros(A) + check_compatible_pattern(A, ag, uplo) if uplo == :F - check_same_pattern(A, S) for k in eachindex(nzA, compressed_indices) nzA[k] = B[compressed_indices[k]] end @@ -534,9 +534,10 @@ end function decompress!( A::AbstractMatrix, B::AbstractMatrix, result::TreeSetColoringResult, uplo::Symbol=:F ) - (; ag, color, reverse_bfs_orders, tree_edge_indices, nt, buffer) = result + (; ag, color, reverse_bfs_orders, tree_edge_indices, nt, diagonal_indices, buffer) = + result (; S) = ag - uplo == :F && check_same_pattern(A, S) + check_compatible_pattern(A, ag, uplo) R = eltype(A) fill!(A, zero(R)) @@ -547,11 +548,9 @@ function decompress!( end # Recover the diagonal coefficients of A - if has_diagonal(ag) - for i in axes(S, 1) - if !iszero(S[i, i]) - A[i, i] = B[i, color[i]] - end + if !augmented_graph(ag) + for i in diagonal_indices + A[i, i] = B[i, color[i]] end end @@ -586,8 +585,23 @@ function decompress!( return A end -function decompress!( - A::SparseMatrixCSC{R}, +""" + decompress_csc!( + nzA::AbstractVector{R}, + A_colptr::AbstractVector, + B::AbstractMatrix{R}, + result::TreeSetColoringResult, + uplo::Symbol=:F, + ) where {R<:Real} + +Decompress the values of `B` into the vector of nonzero entries `nzA` of a +sparse matrix with column pointers `colptr`. This function assumes that the +row indices are sorted in increasing order and are the same as those of the +sparse matrix given to the `coloring` function that returned `result`. +""" +function decompress_csc!( + nzA::AbstractVector{R}, + A_colptr::AbstractVector{<:Integer}, B::AbstractMatrix{R}, result::TreeSetColoringResult, uplo::Symbol=:F, @@ -604,10 +618,6 @@ function decompress!( upper_triangle_offsets, buffer, ) = result - (; S) = ag - A_colptr = A.colptr - nzA = nonzeros(A) - uplo == :F && check_same_pattern(A, S) if eltype(buffer) == R buffer_right_type = buffer @@ -616,7 +626,7 @@ function decompress!( end # Recover the diagonal coefficients of A - if has_diagonal(ag) + if !augmented_graph(ag) if uplo == :L for i in diagonal_indices # A[i, i] is the first element in column i @@ -697,6 +707,17 @@ function decompress!( #! format: on end end + return nothing +end + +function decompress!( + A::SparseMatrixCSC{R}, + B::AbstractMatrix{R}, + result::TreeSetColoringResult, + uplo::Symbol=:F, +) where {R<:Real} + check_compatible_pattern(A, result.ag, uplo) + decompress_csc!(nonzeros(A), A.colptr, B, result, uplo) return A end @@ -708,9 +729,10 @@ function decompress!( result::LinearSystemColoringResult, uplo::Symbol=:F, ) - (; color, strict_upper_nonzero_inds, M_factorization, strict_upper_nonzeros_A) = result - S = result.ag.S - uplo == :F && check_same_pattern(A, S) + (; ag, color, strict_upper_nonzero_inds, M_factorization, strict_upper_nonzeros_A) = + result + S = ag.S + check_compatible_pattern(A, ag, uplo) # TODO: for some reason I cannot use ldiv! with a sparse QR strict_upper_nonzeros_A = M_factorization \ vec(B) @@ -770,10 +792,22 @@ end function decompress!( A::AbstractMatrix, Br::AbstractMatrix, Bc::AbstractMatrix, result::BicoloringResult ) + (; large_colptr, large_rowval, symmetric_result) = result m, n = size(A) + R = eltype(A) + fill!(A, zero(R)) + nzval = Vector{R}(undef, length(large_rowval)) + A_and_noAᵀ = SparseMatrixCSC(m + n, m + n, large_colptr, large_rowval, nzval) Br_and_Bc = _join_compressed!(result, Br, Bc) - A_and_Aᵀ = decompress(Br_and_Bc, result.symmetric_result) - copyto!(A, A_and_Aᵀ[(n + 1):(n + m), 1:n]) # original matrix in bottom left corner + decompress!(A_and_noAᵀ, Br_and_Bc, symmetric_result, :L) + rvA = rowvals(A_and_noAᵀ) + nzA = nonzeros(A_and_noAᵀ) + for j in 1:n + for k in nzrange(A_and_noAᵀ, j) + i = rvA[k] + A[i - n, j] = nzA[k] + end + end return A end @@ -782,10 +816,10 @@ function decompress!( ) (; large_colptr, large_rowval, symmetric_result) = result m, n = size(A) - Br_and_Bc = _join_compressed!(result, Br, Bc) # pretend A is larger A_and_noAᵀ = SparseMatrixCSC(m + n, m + n, large_colptr, large_rowval, A.nzval) # decompress lower triangle only + Br_and_Bc = _join_compressed!(result, Br, Bc) decompress!(A_and_noAᵀ, Br_and_Bc, symmetric_result, :L) return A end diff --git a/src/graph.jl b/src/graph.jl index 9b53a36d..94954602 100644 --- a/src/graph.jl +++ b/src/graph.jl @@ -23,7 +23,7 @@ end SparsityPatternCSC(A::SparseMatrixCSC) = SparsityPatternCSC(A.m, A.n, A.colptr, A.rowval) -Base.eltype(::SparsityPatternCSC{T}) where {T} = T +SparseArrays.indtype(::SparsityPatternCSC{T}) where {T} = T Base.size(S::SparsityPatternCSC) = (S.m, S.n) Base.size(S::SparsityPatternCSC, d::Integer) = d::Integer <= 2 ? size(S)[d] : 1 Base.axes(S::SparsityPatternCSC, d::Integer) = Base.OneTo(size(S, d)) @@ -32,6 +32,11 @@ SparseArrays.nnz(S::SparsityPatternCSC) = length(S.rowval) SparseArrays.rowvals(S::SparsityPatternCSC) = S.rowval SparseArrays.nzrange(S::SparsityPatternCSC, j::Integer) = S.colptr[j]:(S.colptr[j + 1] - 1) +# Needed if using `coloring(::SparsityPatternCSC, ...)` +function Base.similar(A::SparsityPatternCSC, ::Type{T}) where {T} + return SparseArrays.SparseMatrixCSC(A.m, A.n, A.colptr, A.rowval, similar(A.rowval, T)) +end + """ transpose(S::SparsityPatternCSC) @@ -179,6 +184,7 @@ function build_edge_to_index(S::SparsityPatternCSC{T}) where {T} # edge_to_index gives an index for each edge edge_to_index = Vector{T}(undef, nnz(S)) offsets = zeros(T, S.n) + nb_self_loops = 0 counter = 0 rvS = rowvals(S) for j in axes(S, 2) @@ -193,16 +199,17 @@ function build_edge_to_index(S::SparsityPatternCSC{T}) where {T} elseif i == j # this should never be used, make sure it errors edge_to_index[k] = 0 + nb_self_loops += 1 end end end - return edge_to_index + return edge_to_index, nb_self_loops end ## Adjacency graph """ - AdjacencyGraph{T,has_diagonal} + AdjacencyGraph{T,augmented_graph} Undirected graph without self-loops representing the nonzeros of a symmetric matrix (typically a Hessian matrix). @@ -213,45 +220,52 @@ The adjacency graph of a symmetric matrix `A ∈ ℝ^{n × n}` is `G(A) = (V, E) # Constructors - AdjacencyGraph(A::SparseMatrixCSC; has_diagonal::Bool=true) + AdjacencyGraph(A::SparseMatrixCSC; augmented_graph::Bool=false) # Fields -- `S::SparsityPatternCSC{T}`: Underlying sparsity pattern, whose diagonal is empty whenever `has_diagonal` is `false` +- `S::SparsityPatternCSC{T}`: Underlying sparsity pattern, which represents an augmented graph whenever `augmented_graph` is `true`. Here, "augmented graph" means the sparsity pattern of the augmented matrix `H = [0 Jᵀ; J 0]`. - `edge_to_index::Vector{T}`: A vector mapping each nonzero of `S` to a unique edge index (ignoring diagonal and accounting for symmetry, so that `(i, j)` and `(j, i)` get the same index) # References > [_What Color Is Your Jacobian? SparsityPatternCSC Coloring for Computing Derivatives_](https://epubs.siam.org/doi/10.1137/S0036144504444711), Gebremedhin et al. (2005) """ -struct AdjacencyGraph{T<:Integer,has_diagonal} +struct AdjacencyGraph{T<:Integer,augmented_graph} S::SparsityPatternCSC{T} edge_to_index::Vector{T} + nb_self_loops::Int end Base.eltype(::AdjacencyGraph{T}) where {T} = T function AdjacencyGraph( S::SparsityPatternCSC{T}, - edge_to_index::Vector{T}=build_edge_to_index(S); - has_diagonal::Bool=true, + edge_to_index::Vector{T}, + nb_self_loops::Int; + augmented_graph::Bool=false, ) where {T} - return AdjacencyGraph{T,has_diagonal}(S, edge_to_index) + return AdjacencyGraph{T,augmented_graph}(S, edge_to_index, nb_self_loops) +end + +function AdjacencyGraph(S::SparsityPatternCSC; augmented_graph::Bool=false) + edge_to_index, nb_self_loops = build_edge_to_index(S) + return AdjacencyGraph(S, edge_to_index, nb_self_loops; augmented_graph) end -function AdjacencyGraph(A::SparseMatrixCSC; has_diagonal::Bool=true) - return AdjacencyGraph(SparsityPatternCSC(A); has_diagonal) +function AdjacencyGraph(A::SparseMatrixCSC; augmented_graph::Bool=false) + return AdjacencyGraph(SparsityPatternCSC(A); augmented_graph) end -function AdjacencyGraph(A::AbstractMatrix; has_diagonal::Bool=true) - return AdjacencyGraph(SparseMatrixCSC(A); has_diagonal) +function AdjacencyGraph(A::AbstractMatrix; augmented_graph::Bool=false) + return AdjacencyGraph(SparseMatrixCSC(A); augmented_graph) end pattern(g::AdjacencyGraph) = g.S edge_indices(g::AdjacencyGraph) = g.edge_to_index nb_vertices(g::AdjacencyGraph) = pattern(g).n vertices(g::AdjacencyGraph) = 1:nb_vertices(g) -has_diagonal(::AdjacencyGraph{T,hd}) where {T,hd} = hd +augmented_graph(::AdjacencyGraph{T,ag}) where {T,ag} = ag function neighbors(g::AdjacencyGraph, v::Integer) S = pattern(g) @@ -267,30 +281,22 @@ function neighbors_with_edge_indices(g::AdjacencyGraph, v::Integer) return zip(neighbors_v, edges_indices_v) end -degree(g::AdjacencyGraph{T,false}, v::Integer) where {T} = g.S.colptr[v + 1] - g.S.colptr[v] +degree(g::AdjacencyGraph{T,true}, v::Integer) where {T} = g.S.colptr[v + 1] - g.S.colptr[v] -function degree(g::AdjacencyGraph{T,true}, v::Integer) where {T} +function degree(g::AdjacencyGraph{T,false}, v::Integer) where {T} neigh = neighbors(g, v) has_selfloop = insorted(v, neigh) return g.S.colptr[v + 1] - g.S.colptr[v] - has_selfloop end -nb_edges(g::AdjacencyGraph{T,false}) where {T} = nnz(g.S) ÷ 2 - -function nb_edges(g::AdjacencyGraph{T,true}) where {T} - ne = 0 - for v in vertices(g) - ne += degree(g, v) - end - return ne ÷ 2 -end +nb_edges(g::AdjacencyGraph) = (nnz(g.S) - g.nb_self_loops) ÷ 2 maximum_degree(g::AdjacencyGraph) = maximum(Base.Fix1(degree, g), vertices(g)) minimum_degree(g::AdjacencyGraph) = minimum(Base.Fix1(degree, g), vertices(g)) function has_neighbor(g::AdjacencyGraph, v::Integer, u::Integer) for w in neighbors(g, v) - !has_diagonal(g) || (v == w && continue) + augmented_graph(g) || (v == w && continue) if w == u return true end diff --git a/src/interface.jl b/src/interface.jl index 79fe71e9..0d183c9c 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -72,7 +72,7 @@ It is passed as an argument to the main function [`coloring`](@ref). GreedyColoringAlgorithm{decompression}(order=NaturalOrder(); postprocessing=false) GreedyColoringAlgorithm(order=NaturalOrder(); postprocessing=false, decompression=:direct) -- `order::AbstractOrder`: the order in which the columns or rows are colored, which can impact the number of colors. +- `order::Union{AbstractOrder,Tuple}`: the order in which the columns or rows are colored, which can impact the number of colors. Can also be a tuple of different orders to try out, from which the best order (the one with the lowest total number of colors) will be used. - `postprocessing::Bool`: whether or not the coloring will be refined by assigning the neutral color `0` to some vertices. - `decompression::Symbol`: either `:direct` or `:substitution`. Usually `:substitution` leads to fewer colors, at the cost of a more expensive coloring (and decompression). When `:substitution` is not applicable, it falls back on `:direct` decompression. @@ -94,26 +94,31 @@ See their respective docstrings for details. - [`AbstractOrder`](@ref) - [`decompress`](@ref) """ -struct GreedyColoringAlgorithm{decompression,O<:AbstractOrder} <: +struct GreedyColoringAlgorithm{decompression,N,O<:NTuple{N,AbstractOrder}} <: ADTypes.AbstractColoringAlgorithm - order::O + orders::O postprocessing::Bool -end -function GreedyColoringAlgorithm{decompression}( - order::AbstractOrder=NaturalOrder(); postprocessing::Bool=false -) where {decompression} - check_valid_algorithm(decompression) - return GreedyColoringAlgorithm{decompression,typeof(order)}(order, postprocessing) + function GreedyColoringAlgorithm{decompression}( + order_or_orders::Union{AbstractOrder,Tuple}=NaturalOrder(); + postprocessing::Bool=false, + ) where {decompression} + check_valid_algorithm(decompression) + if order_or_orders isa AbstractOrder + orders = (order_or_orders,) + else + orders = order_or_orders + end + return new{decompression,length(orders),typeof(orders)}(orders, postprocessing) + end end function GreedyColoringAlgorithm( - order::AbstractOrder=NaturalOrder(); + order_or_orders::Union{AbstractOrder,Tuple}=NaturalOrder(); postprocessing::Bool=false, decompression::Symbol=:direct, ) - check_valid_algorithm(decompression) - return GreedyColoringAlgorithm{decompression,typeof(order)}(order, postprocessing) + return GreedyColoringAlgorithm{decompression}(order_or_orders; postprocessing) end ## Coloring @@ -225,12 +230,16 @@ function _coloring( ::ColoringProblem{:nonsymmetric,:column}, algo::GreedyColoringAlgorithm, decompression_eltype::Type, - symmetric_pattern::Bool, + symmetric_pattern::Bool; + forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing, ) symmetric_pattern = symmetric_pattern || A isa Union{Symmetric,Hermitian} bg = BipartiteGraph(A; symmetric_pattern) - vertices_in_order = vertices(bg, Val(2), algo.order) - color = partial_distance2_coloring(bg, Val(2), vertices_in_order) + color_by_order = map(algo.orders) do order + vertices_in_order = vertices(bg, Val(2), order) + return partial_distance2_coloring(bg, Val(2), vertices_in_order; forced_colors) + end + color = argmin(maximum, color_by_order) if speed_setting isa WithResult return ColumnColoringResult(A, bg, color) else @@ -244,12 +253,16 @@ function _coloring( ::ColoringProblem{:nonsymmetric,:row}, algo::GreedyColoringAlgorithm, decompression_eltype::Type, - symmetric_pattern::Bool, + symmetric_pattern::Bool; + forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing, ) symmetric_pattern = symmetric_pattern || A isa Union{Symmetric,Hermitian} bg = BipartiteGraph(A; symmetric_pattern) - vertices_in_order = vertices(bg, Val(1), algo.order) - color = partial_distance2_coloring(bg, Val(1), vertices_in_order) + color_by_order = map(algo.orders) do order + vertices_in_order = vertices(bg, Val(1), order) + return partial_distance2_coloring(bg, Val(1), vertices_in_order; forced_colors) + end + color = argmin(maximum, color_by_order) if speed_setting isa WithResult return RowColoringResult(A, bg, color) else @@ -263,11 +276,15 @@ function _coloring( ::ColoringProblem{:symmetric,:column}, algo::GreedyColoringAlgorithm{:direct}, decompression_eltype::Type, - symmetric_pattern::Bool, + symmetric_pattern::Bool; + forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing, ) - ag = AdjacencyGraph(A; has_diagonal=true) - vertices_in_order = vertices(ag, algo.order) - color, star_set = star_coloring(ag, vertices_in_order, algo.postprocessing) + ag = AdjacencyGraph(A; augmented_graph=false) + color_and_star_set_by_order = map(algo.orders) do order + vertices_in_order = vertices(ag, order) + return star_coloring(ag, vertices_in_order, algo.postprocessing; forced_colors) + end + color, star_set = argmin(maximum ∘ first, color_and_star_set_by_order) if speed_setting isa WithResult return StarSetColoringResult(A, ag, color, star_set) else @@ -283,9 +300,18 @@ function _coloring( decompression_eltype::Type{R}, symmetric_pattern::Bool, ) where {R} - ag = AdjacencyGraph(A; has_diagonal=true) - vertices_in_order = vertices(ag, algo.order) - color, tree_set = acyclic_coloring(ag, vertices_in_order, algo.postprocessing) + ag = AdjacencyGraph(A; augmented_graph=false) + color_and_tree_set_by_order = map(algo.orders) do order + vertices_in_order = vertices(ag, order) + return acyclic_coloring(ag, vertices_in_order, algo.postprocessing) + end + # if `color` is empty, `maximum` will fail but `color_and_tree_set_by_order` + # is also one so we can just add a special case for this + if length(color_and_tree_set_by_order) == 1 + color, tree_set = only(color_and_tree_set_by_order) + else + color, tree_set = argmin(maximum ∘ first, color_and_tree_set_by_order) + end if speed_setting isa WithResult return TreeSetColoringResult(A, ag, color, tree_set, R) else @@ -299,19 +325,44 @@ function _coloring( ::ColoringProblem{:nonsymmetric,:bidirectional}, algo::GreedyColoringAlgorithm{:direct}, decompression_eltype::Type{R}, - symmetric_pattern::Bool, + symmetric_pattern::Bool; + forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing, ) where {R} A_and_Aᵀ, edge_to_index = bidirectional_pattern(A; symmetric_pattern) - ag = AdjacencyGraph(A_and_Aᵀ, edge_to_index; has_diagonal=false) - vertices_in_order = vertices(ag, algo.order) - color, star_set = star_coloring(ag, vertices_in_order, algo.postprocessing) + ag = AdjacencyGraph(A_and_Aᵀ, edge_to_index, 0; augmented_graph=true) + outputs_by_order = map(algo.orders) do order + vertices_in_order = vertices(ag, order) + _color, _star_set = star_coloring( + ag, vertices_in_order, algo.postprocessing; forced_colors + ) + (_row_color, _column_color, _symmetric_to_row, _symmetric_to_column) = remap_colors( + eltype(ag), _color, maximum(_color), size(A)... + ) + return ( + _color, + _star_set, + _row_color, + _column_color, + _symmetric_to_row, + _symmetric_to_column, + ) + end + (color, star_set, row_color, column_color, symmetric_to_row, symmetric_to_column) = argmin( + t -> maximum(t[3]) + maximum(t[4]), outputs_by_order + ) # can't use ncolors without computing the full result if speed_setting isa WithResult symmetric_result = StarSetColoringResult(A_and_Aᵀ, ag, color, star_set) - return BicoloringResult(A, ag, symmetric_result, R) - else - row_color, column_color, _ = remap_colors( - eltype(ag), color, maximum(color), size(A)... + return BicoloringResult( + A, + ag, + symmetric_result, + row_color, + column_color, + symmetric_to_row, + symmetric_to_column, + R, ) + else return row_color, column_color end end @@ -325,16 +376,38 @@ function _coloring( symmetric_pattern::Bool, ) where {R} A_and_Aᵀ, edge_to_index = bidirectional_pattern(A; symmetric_pattern) - ag = AdjacencyGraph(A_and_Aᵀ, edge_to_index; has_diagonal=false) - vertices_in_order = vertices(ag, algo.order) - color, tree_set = acyclic_coloring(ag, vertices_in_order, algo.postprocessing) + ag = AdjacencyGraph(A_and_Aᵀ, edge_to_index, 0; augmented_graph=true) + outputs_by_order = map(algo.orders) do order + vertices_in_order = vertices(ag, order) + _color, _tree_set = acyclic_coloring(ag, vertices_in_order, algo.postprocessing) + (_row_color, _column_color, _symmetric_to_row, _symmetric_to_column) = remap_colors( + eltype(ag), _color, maximum(_color), size(A)... + ) + return ( + _color, + _tree_set, + _row_color, + _column_color, + _symmetric_to_row, + _symmetric_to_column, + ) + end + (color, tree_set, row_color, column_color, symmetric_to_row, symmetric_to_column) = argmin( + t -> maximum(t[3]) + maximum(t[4]), outputs_by_order + ) # can't use ncolors without computing the full result if speed_setting isa WithResult symmetric_result = TreeSetColoringResult(A_and_Aᵀ, ag, color, tree_set, R) - return BicoloringResult(A, ag, symmetric_result, R) - else - row_color, column_color, _ = remap_colors( - eltype(ag), color, maximum(color), size(A)... + return BicoloringResult( + A, + ag, + symmetric_result, + row_color, + column_color, + symmetric_to_row, + symmetric_to_column, + R, ) + else return row_color, column_color end end diff --git a/src/matrices.jl b/src/matrices.jl index eaf79882..453e3062 100644 --- a/src/matrices.jl +++ b/src/matrices.jl @@ -62,24 +62,57 @@ function respectful_similar(A::Union{Symmetric,Hermitian}, ::Type{T}) where {T} end """ - same_pattern(A, B) + compatible_pattern(A::AbstractMatrix, bg::BipartiteGraph) + compatible_pattern(A::AbstractMatrix, ag::AdjacencyGraph, uplo::Symbol) -Perform a partial equality check on the sparsity patterns of `A` and `B`: +Perform a coarse compatibility check between the sparsity pattern of `A` +and the reference sparsity pattern encoded in a graph structure. -- if the return is `true`, they might have the same sparsity pattern but we're not sure -- if the return is `false`, they definitely don't have the same sparsity pattern +This function only checks necessary conditions for the two sparsity patterns to match. +In particular, it may return `true` even if the patterns are not identical. + +When A is a `SparseMatrixCSC`, additional checks on the sparsity structure are performed. + +# Return value +- `true` : the sparsity patterns are potentially compatible +- `false` : the sparsity patterns are definitely incompatible """ -same_pattern(A, B) = size(A) == size(B) +compatible_pattern(A::AbstractMatrix, bg::BipartiteGraph) = size(A) == size(bg.S2) +function compatible_pattern(A::SparseMatrixCSC, bg::BipartiteGraph) + return size(A) == size(bg.S2) && nnz(A) == nnz(bg.S2) +end + +function compatible_pattern(A::AbstractMatrix, ag::AdjacencyGraph, uplo::Symbol) + return size(A) == size(ag.S) +end + +function compatible_pattern(A::SparseMatrixCSC, ag::AdjacencyGraph, uplo::Symbol) + nnzS = (uplo == :L || uplo == :U) ? (nb_edges(ag) + ag.nb_self_loops) : nnz(ag.S) + return size(A) == size(ag.S) && nnz(A) == nnzS +end -function same_pattern( - A::Union{SparseMatrixCSC,SparsityPatternCSC}, - B::Union{SparseMatrixCSC,SparsityPatternCSC}, -) - return size(A) == size(B) && nnz(A) == nnz(B) +function check_compatible_pattern(A::AbstractMatrix, bg::BipartiteGraph) + if !compatible_pattern(A, bg) + throw(DimensionMismatch("`A` and `bg.S2` must have the same sparsity pattern.")) + end end -function check_same_pattern(A, S) - if !same_pattern(A, S) - throw(DimensionMismatch("`A` and `S` must have the same sparsity pattern.")) +function check_compatible_pattern(A::AbstractMatrix, ag::AdjacencyGraph, uplo::Symbol) + if !compatible_pattern(A, ag, uplo) + if uplo == :L + throw( + DimensionMismatch( + "`A` and `tril(ag.S)` must have the same sparsity pattern." + ), + ) + elseif uplo == :U + throw( + DimensionMismatch( + "`A` and `triu(ag.S)` must have the same sparsity pattern." + ), + ) + else # uplo == :F + throw(DimensionMismatch("`A` and `ag.S` must have the same sparsity pattern.")) + end end end diff --git a/src/optimal.jl b/src/optimal.jl new file mode 100644 index 00000000..2a70a59d --- /dev/null +++ b/src/optimal.jl @@ -0,0 +1,34 @@ +""" + OptimalColoringAlgorithm + +Coloring algorithm that relies on mathematical programming with [JuMP](https://jump.dev/) to find an optimal coloring. + +!!! warning + This algorithm is only available when JuMP is loaded. If you encounter a method error, run `import JuMP` in your REPL and try again. + It only works for nonsymmetric, unidirectional colorings problems. + +!!! danger + The coloring problem is NP-hard, so it is unreasonable to expect an optimal solution in reasonable time for large instances. + +# Constructor + + OptimalColoringAlgorithm(optimizer; silent::Bool=true, assert_solved::Bool=true) + +The `optimizer` argument can be any JuMP-compatible optimizer. +However, the problem formulation is best suited to CP-SAT optimizers like [MiniZinc](https://github.com/jump-dev/MiniZinc.jl). +You can use [`optimizer_with_attributes`](https://jump.dev/JuMP.jl/stable/api/JuMP/#optimizer_with_attributes) to set solver-specific parameters. + +# Keyword arguments + +- `silent`: whether to suppress solver output +- `assert_solved`: whether to check that the solver found an optimal solution (as opposed to running out of time for example) +""" +struct OptimalColoringAlgorithm{O} <: ADTypes.AbstractColoringAlgorithm + optimizer::O + silent::Bool + assert_solved::Bool +end + +function OptimalColoringAlgorithm(optimizer; silent::Bool=true, assert_solved::Bool=true) + return OptimalColoringAlgorithm(optimizer, silent, assert_solved) +end diff --git a/src/order.jl b/src/order.jl index aab5bb23..20884c6a 100644 --- a/src/order.jl +++ b/src/order.jl @@ -334,7 +334,7 @@ function vertices( π[index] = u for v in neighbors(g, u) - !has_diagonal(g) || (u == v && continue) + augmented_graph(g) || (u == v && continue) dv = degrees[v] dv == -1 && continue update_bucket!(db, v, dv; degtype, direction) diff --git a/src/postprocessing.jl b/src/postprocessing.jl new file mode 100644 index 00000000..7b2b58f5 --- /dev/null +++ b/src/postprocessing.jl @@ -0,0 +1,168 @@ +## Postprocessing + +function postprocess!( + color::AbstractVector{<:Integer}, + star_or_tree_set::Union{StarSet,TreeSet}, + g::AdjacencyGraph, + offsets::AbstractVector{<:Integer}, +) + S = pattern(g) + edge_to_index = edge_indices(g) + # flag which colors are actually used during decompression + nb_colors = maximum(color) + color_used = zeros(Bool, nb_colors) + + # nonzero diagonal coefficients force the use of their respective color (there can be no neutral colors if the diagonal is fully nonzero) + if !augmented_graph(g) + for i in axes(S, 1) + if !iszero(S[i, i]) + color_used[color[i]] = true + end + end + end + + if star_or_tree_set isa StarSet + # star_or_tree_set is a StarSet + postprocess_with_star_set!(g, color_used, color, star_or_tree_set) + else + # star_or_tree_set is a TreeSet + postprocess_with_tree_set!(color_used, color, star_or_tree_set) + end + + # if at least one of the colors is useless, modify the color assignments of vertices + if any(!, color_used) + num_colors_useless = 0 + + # determine what are the useless colors and compute the offsets + for ci in 1:nb_colors + if color_used[ci] + offsets[ci] = num_colors_useless + else + num_colors_useless += 1 + end + end + + # assign the neutral color to every vertex with a useless color and remap the colors + for i in eachindex(color) + ci = color[i] + if !color_used[ci] + # assign the neutral color + color[i] = 0 + else + # remap the color to not have any gap + color[i] -= offsets[ci] + end + end + end + return color +end + +function postprocess_with_star_set!( + g::AdjacencyGraph, + color_used::Vector{Bool}, + color::AbstractVector{<:Integer}, + star_set::StarSet, +) + S = pattern(g) + edge_to_index = edge_indices(g) + + # only the colors of the hubs are used + (; star, hub) = star_set + nb_trivial_stars = 0 + + # Iterate through all non-trivial stars + for s in eachindex(hub) + h = hub[s] + if h > 0 + color_used[color[h]] = true + else + nb_trivial_stars += 1 + end + end + + # Process the trivial stars (if any) + if nb_trivial_stars > 0 + rvS = rowvals(S) + for j in axes(S, 2) + for k in nzrange(S, j) + i = rvS[k] + if i > j + index_ij = edge_to_index[k] + s = star[index_ij] + h = hub[s] + if h < 0 + h = abs(h) + spoke = h == j ? i : j + if color_used[color[spoke]] + # Switch the hub and the spoke to possibly avoid adding one more used color + hub[s] = spoke + else + # Keep the current hub + color_used[color[h]] = true + end + end + end + end + end + end + return color_used +end + +function postprocess_with_tree_set!( + color_used::Vector{Bool}, color::AbstractVector{<:Integer}, tree_set::TreeSet +) + # only the colors of non-leaf vertices are used + (; reverse_bfs_orders, is_star, tree_edge_indices, nt) = tree_set + nb_trivial_trees = 0 + + # Iterate through all non-trivial trees + for k in 1:nt + # Position of the first edge in the tree + first = tree_edge_indices[k] + + # Total number of edges in the tree + ne_tree = tree_edge_indices[k + 1] - first + + # Check if we have more than one edge in the tree (non-trivial tree) + if ne_tree > 1 + # Determine if the tree is a star + if is_star[k] + # It is a non-trivial star and only the color of the hub is needed + (_, hub) = reverse_bfs_orders[first] + color_used[color[hub]] = true + else + # It is not a star and both colors are needed during the decompression + (i, j) = reverse_bfs_orders[first] + color_used[color[i]] = true + color_used[color[j]] = true + end + else + nb_trivial_trees += 1 + end + end + + # Process the trivial trees (if any) + if nb_trivial_trees > 0 + for k in 1:nt + # Position of the first edge in the tree + first = tree_edge_indices[k] + + # Total number of edges in the tree + ne_tree = tree_edge_indices[k + 1] - first + + # Check if we have exactly one edge in the tree + if ne_tree == 1 + (i, j) = reverse_bfs_orders[first] + if color_used[color[i]] + # Make i the root to avoid possibly adding one more used color + # Switch it with the (only) leaf + reverse_bfs_orders[first] = (j, i) + else + # Keep j as the root + color_used[color[j]] = true + end + end + end + end + return color_used +end diff --git a/src/result.jl b/src/result.jl index b5b9cdd2..0d47cf13 100644 --- a/src/result.jl +++ b/src/result.jl @@ -82,6 +82,9 @@ Create a color-indexed vector `group` such that `i ∈ group[c]` iff `color[i] = Assumes the colors are contiguously numbered from `0` to some `cmax`. """ function group_by_color(::Type{T}, color::AbstractVector) where {T<:Integer} + if isempty(color) + return typeof(view(T[], 1:0))[] + end cmin, cmax = extrema(color) @assert cmin >= 0 # Compute group sizes and offsets for a joint storage @@ -401,31 +404,31 @@ function TreeSetColoringResult( decompression_eltype::Type{R}, ) where {T<:Integer,R} (; reverse_bfs_orders, tree_edge_indices, nt) = tree_set - (; S) = ag + (; S, nb_self_loops) = ag nvertices = length(color) group = group_by_color(T, color) rv = rowvals(S) # Vector for the decompression of the diagonal coefficients - diagonal_indices = T[] - diagonal_nzind = T[] - ndiag = 0 + diagonal_indices = Vector{T}(undef, nb_self_loops) + diagonal_nzind = Vector{T}(undef, nb_self_loops) - if has_diagonal(ag) + if !augmented_graph(ag) + l = 0 for j in axes(S, 2) for k in nzrange(S, j) i = rv[k] if i == j - push!(diagonal_indices, i) - push!(diagonal_nzind, k) - ndiag += 1 + l += 1 + diagonal_indices[l] = i + diagonal_nzind[l] = k end end end end # Vectors for the decompression of the off-diagonal coefficients - nedges = (nnz(S) - ndiag) ÷ 2 + nedges = nb_edges(ag) lower_triangle_offsets = Vector{T}(undef, nedges) upper_triangle_offsets = Vector{T}(undef, nedges) @@ -686,14 +689,15 @@ function BicoloringResult( A::AbstractMatrix, ag::AdjacencyGraph{T}, symmetric_result::AbstractColoringResult{:symmetric,:column}, + row_color::Vector{T}, + column_color::Vector{T}, + symmetric_to_row::Vector{T}, + symmetric_to_column::Vector{T}, decompression_eltype::Type{R}, ) where {T,R} m, n = size(A) symmetric_color = column_colors(symmetric_result) num_sym_colors = maximum(symmetric_color) - row_color, column_color, symmetric_to_row, symmetric_to_column = remap_colors( - T, symmetric_color, num_sym_colors, m, n - ) column_group = group_by_color(T, column_color) row_group = group_by_color(T, row_color) Br_and_Bc = Matrix{R}(undef, n + m, num_sym_colors) diff --git a/src/structured.jl b/src/structured.jl new file mode 100644 index 00000000..812256bf --- /dev/null +++ b/src/structured.jl @@ -0,0 +1,181 @@ +""" + StructuredColoringAlgorithm <: ADTypes.AbstractColoringAlgorithm + +Coloring algorithm which leverages specific matrix structures to produce optimal or near-optimal solutions. + +The following matrix types are supported: + +- From the standard library `LinearAlgebra`: `Diagonal`, `Bidiagonal`, `Tridiagonal` +- From [BandedMatrices.jl](https://github.com/JuliaLinearAlgebra/BandedMatrices.jl): [`BandedMatrix`](@extref BandedMatrices.BandedMatrix) + +!!! warning + Only `:nonsymmetric` structures with `:row` or `:column` partitions (aka unidirectional Jacobian colorings) are supported by this algorithm at the moment. + +!!! tip + To request support for a new type of structured matrix, open an issue on the SparseMatrixColorings.jl GitHub repository! +""" +struct StructuredColoringAlgorithm <: ADTypes.AbstractColoringAlgorithm end + +#= +This code is partly taken from ArrayInterface.jl +https://github.com/JuliaArrays/ArrayInterface.jl +=# + +""" + cycle_range(k::Integer, n::Integer) + +Concatenate copies of `1:k` to fill a vector of length `n` (with one partial copy allowed at the end). +""" +function cycle_range(k::Integer, n::Integer) + color = Vector{Int}(undef, n) + for i in eachindex(color) + color[i] = 1 + (i - 1) % k + end + return color +end + +## Diagonal + +function coloring( + A::Diagonal, + ::ColoringProblem{:nonsymmetric,:column}, + ::StructuredColoringAlgorithm; + kwargs..., +) + color = fill(1, size(A, 2)) + bg = BipartiteGraph(A) + return ColumnColoringResult(A, bg, color) +end + +function coloring( + A::Diagonal, + ::ColoringProblem{:nonsymmetric,:row}, + ::StructuredColoringAlgorithm; + kwargs..., +) + color = fill(1, size(A, 1)) + bg = BipartiteGraph(A) + return RowColoringResult(A, bg, color) +end + +function decompress!(A::Diagonal, B::AbstractMatrix, result::ColumnColoringResult) + color = column_colors(result) + for j in axes(A, 2) + A[j, j] = B[j, color[j]] + end + return A +end + +function decompress!(A::Diagonal, B::AbstractMatrix, result::RowColoringResult) + color = row_colors(result) + for i in axes(A, 1) + A[i, i] = B[color[i], i] + end + return A +end + +## Bidiagonal + +function coloring( + A::Bidiagonal, + ::ColoringProblem{:nonsymmetric,:column}, + ::StructuredColoringAlgorithm; + kwargs..., +) + color = cycle_range(2, size(A, 2)) + bg = BipartiteGraph(A) + return ColumnColoringResult(A, bg, color) +end + +function coloring( + A::Bidiagonal, + ::ColoringProblem{:nonsymmetric,:row}, + ::StructuredColoringAlgorithm; + kwargs..., +) + color = cycle_range(2, size(A, 1)) + bg = BipartiteGraph(A) + return RowColoringResult(A, bg, color) +end + +function decompress!(A::Bidiagonal, B::AbstractMatrix, result::ColumnColoringResult) + color = column_colors(result) + for j in axes(A, 2) + c = color[j] + A[j, j] = B[j, c] + if A.uplo == 'U' && j > 1 # above + A[j - 1, j] = B[j - 1, c] + elseif A.uplo == 'L' && j < size(A, 2) # below + A[j + 1, j] = B[j + 1, c] + end + end + return A +end + +function decompress!(A::Bidiagonal, B::AbstractMatrix, result::RowColoringResult) + color = row_colors(result) + for i in axes(A, 1) + c = color[i] + A[i, i] = B[c, i] + if A.uplo == 'U' && i < size(A, 1) # right + A[i, i + 1] = B[c, i + 1] + elseif A.uplo == 'L' && i > 1 # left + A[i, i - 1] = B[c, i - 1] + end + end + return A +end + +## Tridiagonal + +function coloring( + A::Tridiagonal, + ::ColoringProblem{:nonsymmetric,:column}, + ::StructuredColoringAlgorithm; + kwargs..., +) + color = cycle_range(3, size(A, 2)) + bg = BipartiteGraph(A) + return ColumnColoringResult(A, bg, color) +end + +function coloring( + A::Tridiagonal, + ::ColoringProblem{:nonsymmetric,:row}, + ::StructuredColoringAlgorithm; + kwargs..., +) + color = cycle_range(3, size(A, 1)) + bg = BipartiteGraph(A) + return RowColoringResult(A, bg, color) +end + +function decompress!(A::Tridiagonal, B::AbstractMatrix, result::ColumnColoringResult) + color = column_colors(result) + for j in axes(A, 2) + c = color[j] + A[j, j] = B[j, c] + if j > 1 # above + A[j - 1, j] = B[j - 1, c] + end + if j < size(A, 2) # below + A[j + 1, j] = B[j + 1, c] + end + end + return A +end + +function decompress!(A::Tridiagonal, B::AbstractMatrix, result::RowColoringResult) + color = row_colors(result) + for i in axes(A, 1) + c = color[i] + A[i, i] = B[c, i] + if i < size(A, 1) # right + A[i, i + 1] = B[c, i + 1] + end + if i > 1 # left + A[i, i - 1] = B[c, i - 1] + end + end + return A +end diff --git a/test/Project.toml b/test/Project.toml index ab595aba..5e41ceb2 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -12,12 +12,15 @@ CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8" Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" -JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899" +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MatrixDepot = "b51810bb-c9f3-55da-ae3c-350fc1fbce05" +MiniZinc = "a7f392d2-6c35-496e-b8cc-0974fbfcbf91" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +cuSPARSE = "b26da814-b3bc-49ef-b0ee-c816305aa060" diff --git a/test/adtypes.jl b/test/adtypes.jl index bbb124a6..da0807d3 100644 --- a/test/adtypes.jl +++ b/test/adtypes.jl @@ -1,26 +1,49 @@ using ADTypes: ADTypes using SparseArrays +using LinearAlgebra using SparseMatrixColorings using Test -@testset "Column coloring" begin - problem = ColoringProblem(; structure=:nonsymmetric, partition=:column) - algo = ADTypes.NoColoringAlgorithm() - A = sprand(10, 20, 0.1) - result = coloring(A, problem, algo) - B = compress(A, result) - @test size(B) == size(A) - @test column_colors(result) == ADTypes.column_coloring(A, algo) - @test decompress(B, result) == A -end +@testset "NoColoringAlgorithm" begin + @testset "Column coloring" begin + problem = ColoringProblem(; structure=:nonsymmetric, partition=:column) + algo = ADTypes.NoColoringAlgorithm() + A = sprand(10, 20, 0.3) + result = coloring(A, problem, algo) + B = compress(A, result) + @test size(B) == size(A) + @test column_colors(result) == ADTypes.column_coloring(A, algo) + @test decompress(B, result) == A + end + + @testset "Row coloring" begin + problem = ColoringProblem(; structure=:nonsymmetric, partition=:row) + algo = ADTypes.NoColoringAlgorithm() + A = sprand(10, 20, 0.3) + result = coloring(A, problem, algo) + B = compress(A, result) + @test size(B) == size(A) + @test row_colors(result) == ADTypes.row_coloring(A, algo) + @test decompress(B, result) == A + end + + @testset "Symmetric coloring" begin + problem = ColoringProblem(; structure=:symmetric, partition=:column) + algo = ADTypes.NoColoringAlgorithm() + A = Symmetric(sprand(20, 20, 0.3)) + result = coloring(A, problem, algo) + B = compress(A, result) + @test size(B) == size(A) + @test column_colors(result) == ADTypes.column_coloring(A, algo) + @test decompress(B, result) == A + end -@testset "Row coloring" begin - problem = ColoringProblem(; structure=:nonsymmetric, partition=:row) - algo = ADTypes.NoColoringAlgorithm() - A = sprand(10, 20, 0.1) - result = coloring(A, problem, algo) - B = compress(A, result) - @test size(B) == size(A) - @test row_colors(result) == ADTypes.row_coloring(A, algo) - @test decompress(B, result) == A + @testset "Bicoloring" begin + problem = ColoringProblem(; structure=:nonsymmetric, partition=:bidirectional) + algo = ADTypes.NoColoringAlgorithm() + A = sprand(10, 20, 0.3) + result = coloring(A, problem, algo) + Br, Bc = compress(A, result) + @test decompress(Br, Bc, result) == A + end end diff --git a/test/allocations.jl b/test/allocations.jl index caed10f8..31305e63 100644 --- a/test/allocations.jl +++ b/test/allocations.jl @@ -154,7 +154,7 @@ end b64 = @b fast_coloring(A64, problem, algo) b32 = @b fast_coloring(A32, problem, algo) # check that we allocate no more than 50% + epsilon with Int32 - @test b32.bytes < 0.6 * b64.bytes + @test b32.bytes < 0.7 * b64.bytes end end end; diff --git a/test/check.jl b/test/check.jl index 4e9e5a61..8cf2c495 100644 --- a/test/check.jl +++ b/test/check.jl @@ -1,9 +1,12 @@ using LinearAlgebra +using SparseArrays using SparseMatrixColorings: structurally_orthogonal_columns, symmetrically_orthogonal_columns, structurally_biorthogonal, directly_recoverable_columns, + substitutable_columns, + substitutable_bidirectional, what_fig_41, efficient_fig_1 using Test @@ -184,3 +187,135 @@ For coefficient (i=1, j=1) with colors (ci=0, cj=0): """, ) structurally_biorthogonal(A, [0, 2, 2, 3], [0, 2, 2, 2, 3], verbose=true) end + +@testset "Substitutable columns" begin + A1 = [ + 1 1 1 1 1 + 1 1 0 0 0 + 1 0 1 0 0 + 1 0 0 1 0 + 1 0 0 0 1 + ] + B1 = [ + 1 6 7 8 9 + 6 2 0 0 0 + 7 0 3 0 0 + 8 0 0 4 0 + 9 0 0 0 5 + ] + A2 = [ + 1 1 0 0 0 + 1 1 1 0 0 + 0 1 1 1 0 + 0 0 1 1 1 + 0 0 0 1 1 + ] + B2 = [ + 5 1 0 0 0 + 1 6 2 0 0 + 0 2 7 3 0 + 0 0 3 8 4 + 0 0 0 4 9 + ] + A3 = [ + 0 1 1 1 1 + 1 0 1 1 1 + 1 1 0 1 1 + 1 1 1 0 1 + 1 1 1 1 0 + ] + B3 = [ + 0 1 2 3 4 + 1 0 5 6 7 + 2 5 0 8 9 + 3 6 8 0 10 + 4 7 9 10 0 + ] + + # success + + @test substitutable_columns(A1, B1, [1, 2, 2, 2, 2]) + @test substitutable_columns(A2, B2, [1, 2, 3, 1, 2]) + @test substitutable_columns(A3, B3, [1, 2, 3, 4, 0]) + + # failure + + @test !substitutable_columns(A1, B1, [1, 1, 1, 1]) + log = (:warn, "4 colors provided for 5 columns.") + @test_logs log substitutable_columns(A1, B1, [1, 1, 1, 1]; verbose=true) + + @test !substitutable_columns(A1, B1, [1, 1, 1, 1, 1]) + @test_logs ( + :warn, + """ +For coefficient (i=1, j=1) with colors (ci=1, cj=1): +- For the row 5 in row color ci=1, A[5, 1] is ordered after A[1, 1]. +- For the column 5 in column color cj=1, A[1, 5] is ordered after A[1, 1]. +""", + ) substitutable_columns(A1, B1, [1, 1, 1, 1, 1]; verbose=true) + + @test !substitutable_columns(A2, B2, [1, 2, 0, 1, 2]) + @test_logs ( + :warn, + """ +For coefficient (i=3, j=3) with colors (ci=0, cj=0): +- Row color ci=0 is neutral. +- Column color cj=0 is neutral. +""", + ) substitutable_columns(A2, B2, [1, 2, 0, 1, 2]; verbose=true) + + @test !substitutable_columns(A3, B3, [0, 1, 2, 3, 3]) + @test_logs ( + :warn, + """ +For coefficient (i=1, j=4) with colors (ci=0, cj=3): +- Row color ci=0 is neutral. +- For the column 5 in column color cj=3, A[1, 5] is ordered after A[1, 4]. +""", + ) substitutable_columns(A3, B3, [0, 1, 2, 3, 3]; verbose=true) + + @test !substitutable_columns(A3, B3, [1, 2, 3, 3, 0]) + @test_logs ( + :warn, + """ +For coefficient (i=3, j=5) with colors (ci=3, cj=0): +- For the row 4 in row color ci=3, A[4, 5] is ordered after A[3, 5]. +- Column color cj=0 is neutral. +""", + ) substitutable_columns(A3, B3, [1, 2, 3, 3, 0]; verbose=true) +end + +@testset "Substitutable bidirectional" begin + A = [ + 1 0 0 + 0 1 0 + 0 0 1 + ] + B = [ + 1 0 0 + 0 2 0 + 0 0 3 + ] + + # success + + @test substitutable_bidirectional(A, B, [1, 0, 0], [0, 1, 1]) + + # failure + + log = (:warn, "2 colors provided for 3 columns.") + @test_logs log !substitutable_bidirectional(A, B, [1, 0, 0], [0, 1]; verbose=true) + + log = (:warn, "4 colors provided for 3 rows.") + @test_logs log !substitutable_bidirectional(A, B, [1, 0, 0, 1], [0, 1, 1]; verbose=true) +end + +# See https://github.com/JuliaDiff/SparseMatrixColorings.jl/pull/300 +@testset "Empty matrix" begin + problem = ColoringProblem(; structure=:symmetric, partition=:column) + algo = GreedyColoringAlgorithm(; decompression=:substitution) + S = spzeros(Int, 0, 0) + result = coloring(S, problem, algo) + @test isempty(result.color) + @test isempty(result.group) +end diff --git a/test/constant.jl b/test/constant.jl index 5de881b8..97c7f238 100644 --- a/test/constant.jl +++ b/test/constant.jl @@ -1,16 +1,20 @@ using ADTypes: ADTypes using SparseMatrixColorings +using SparseMatrixColorings: InvalidColoringError using Test -matrix_template = ones(100, 200) +matrix_template = ones(Bool, 10, 20) +sym_matrix_template = ones(Bool, 10, 10) @testset "Column coloring" begin problem = ColoringProblem(; structure=:nonsymmetric, partition=:column) - color = rand(1:5, size(matrix_template, 2)) + color = collect(1:20) algo = ConstantColoringAlgorithm(matrix_template, color; partition=:column) - wrong_algo = ConstantColoringAlgorithm(matrix_template, color; partition=:row) + wrong_algo = ConstantColoringAlgorithm{:row}(matrix_template, color) + wrong_color = ConstantColoringAlgorithm{:column}(matrix_template, ones(Int, 20)) @test_throws DimensionMismatch coloring(transpose(matrix_template), problem, algo) @test_throws MethodError coloring(matrix_template, problem, wrong_algo) + @test_throws InvalidColoringError coloring(matrix_template, problem, wrong_color) result = coloring(matrix_template, problem, algo) @test column_colors(result) == color @test ADTypes.column_coloring(matrix_template, algo) == color @@ -19,9 +23,13 @@ end @testset "Row coloring" begin problem = ColoringProblem(; structure=:nonsymmetric, partition=:row) - color = rand(1:5, size(matrix_template, 1)) + color = collect(1:10) algo = ConstantColoringAlgorithm(matrix_template, color; partition=:row) + wrong_algo = ConstantColoringAlgorithm{:column}(matrix_template, color) + wrong_color = ConstantColoringAlgorithm{:row}(matrix_template, ones(Int, 10)) @test_throws DimensionMismatch coloring(transpose(matrix_template), problem, algo) + @test_throws MethodError coloring(matrix_template, problem, wrong_algo) + @test_throws InvalidColoringError coloring(matrix_template, problem, wrong_color) result = coloring(matrix_template, problem, algo) @test row_colors(result) == color @test ADTypes.row_coloring(matrix_template, algo) == color @@ -29,9 +37,22 @@ end end @testset "Symmetric coloring" begin - wrong_problem = ColoringProblem(; structure=:symmetric, partition=:column) - color = rand(1:5, size(matrix_template, 2)) - algo = ConstantColoringAlgorithm(matrix_template, color; partition=:column) - @test_throws MethodError coloring(matrix_template, wrong_problem, algo) - @test_throws MethodError ADTypes.symmetric_coloring(matrix_template, algo) + problem = ColoringProblem(; structure=:symmetric, partition=:column) + color = collect(1:10) + algo = ConstantColoringAlgorithm( + sym_matrix_template, color; partition=:column, structure=:symmetric + ) + wrong_algo = ConstantColoringAlgorithm{:column,:nonsymmetric}( + sym_matrix_template, color + ) + wrong_color = ConstantColoringAlgorithm{:column,:symmetric}( + sym_matrix_template, ones(Int, 20) + ) + @test_throws DimensionMismatch coloring(matrix_template, problem, algo) + @test_throws MethodError coloring(sym_matrix_template, problem, wrong_algo) + @test_throws InvalidColoringError coloring(sym_matrix_template, problem, wrong_color) + result = coloring(sym_matrix_template, problem, algo) + @test column_colors(result) == color + @test ADTypes.symmetric_coloring(sym_matrix_template, algo) == color + @test_throws MethodError ADTypes.column_coloring(sym_matrix_template, algo) end diff --git a/test/cuda.jl b/test/cuda.jl index 6549b168..8b27e07b 100644 --- a/test/cuda.jl +++ b/test/cuda.jl @@ -1,4 +1,4 @@ -using CUDA.CUSPARSE: CuSparseMatrixCSC, CuSparseMatrixCSR +using cuSPARSE: CuSparseMatrixCSC, CuSparseMatrixCSR using LinearAlgebra using SparseArrays using SparseMatrixColorings @@ -6,8 +6,6 @@ import SparseMatrixColorings as SMC using StableRNGs using Test -include("utils.jl") - rng = StableRNG(63) asymmetric_params = vcat( diff --git a/test/graph.jl b/test/graph.jl index 9c8f66c1..773a0206 100644 --- a/test/graph.jl +++ b/test/graph.jl @@ -15,7 +15,8 @@ using Test ## SparsityPatternCSC @testset "SparsityPatternCSC" begin - @test eltype(SparsityPatternCSC(sprand(10, 10, 0.1))) == Int + @test eltype(SparsityPatternCSC(sprand(10, 10, 0.1))) == Bool + @test SparseArrays.indtype(SparsityPatternCSC(sprand(10, 10, 0.1))) == Int @testset "Transpose" begin for _ in 1:1000 m, n = rand(100:1000), rand(100:1000) @@ -166,7 +167,7 @@ end; @test degree(g, 7) == 6 @test degree(g, 8) == 6 - g = AdjacencyGraph(transpose(A) * A; has_diagonal=false) + g = AdjacencyGraph(transpose(A) * A; augmented_graph=true) # wrong degree @test degree(g, 1) == 4 @test degree(g, 2) == 4 diff --git a/test/matrices.jl b/test/matrices.jl index 1571497b..fa76946c 100644 --- a/test/matrices.jl +++ b/test/matrices.jl @@ -1,7 +1,12 @@ using LinearAlgebra using SparseArrays using SparseMatrixColorings: - check_same_pattern, matrix_versions, respectful_similar, same_pattern + BipartiteGraph, + AdjacencyGraph, + matrix_versions, + respectful_similar, + compatible_pattern, + check_compatible_pattern using StableRNGs using Test @@ -33,7 +38,7 @@ same_view(::Adjoint, ::Adjoint) = true end end -@testset "Sparsity pattern" begin +@testset "Compatible sparsity pattern -- BipartiteGraph" begin S = sparse([ 0 1 1 0 1 0 @@ -44,9 +49,33 @@ end A2 = copy(S) A2[1, 1] = 1 - @test same_pattern(A1, S) - @test !same_pattern(A2, S) - @test same_pattern(Matrix(A2), S) + bg1 = BipartiteGraph(A1) + bg2 = BipartiteGraph(A2) + @test compatible_pattern(S, bg1) + @test !compatible_pattern(S, bg2) + @test compatible_pattern(Matrix(S), bg1) - @test_throws DimensionMismatch check_same_pattern(A2, S) + @test_throws DimensionMismatch check_compatible_pattern(S, bg2) +end + +@testset "Compatible sparsity pattern -- AdjacencyGraph" begin + S = sparse([ + 1 0 1 + 0 1 1 + 1 1 0 + ]) + + A1 = copy(S) + A2 = copy(S) + A2[3, 3] = 1 + + ag1 = AdjacencyGraph(A1) + ag2 = AdjacencyGraph(A2) + for (op, uplo) in ((tril, :L), (triu, :U), (identity, :F)) + @test compatible_pattern(op(S), ag1, uplo) + @test !compatible_pattern(op(S), ag2, uplo) + @test compatible_pattern(Matrix(op(S)), ag1, uplo) + + @test_throws DimensionMismatch check_compatible_pattern(op(S), ag2, uplo) + end end diff --git a/test/optimal.jl b/test/optimal.jl new file mode 100644 index 00000000..36bf91c6 --- /dev/null +++ b/test/optimal.jl @@ -0,0 +1,46 @@ +using SparseArrays +using SparseMatrixColorings +using StableRNGs +using Test +using JuMP +using MiniZinc +using HiGHS + +rng = StableRNG(0) + +asymmetric_params = vcat( + [(10, 20, p) for p in (0.0:0.1:0.5)], [(20, 10, p) for p in (0.0:0.1:0.5)] +) + +algo = GreedyColoringAlgorithm() +optalgo = OptimalColoringAlgorithm(() -> MiniZinc.Optimizer{Float64}("highs"); silent=false) + +# TODO: reactivate tests once https://github.com/jump-dev/MiniZinc.jl/issues/103 is fixed + +@testset "Column coloring" begin + problem = ColoringProblem(; structure=:nonsymmetric, partition=:column) + for (m, n, p) in asymmetric_params + A = sprand(rng, m, n, p) + result = coloring(A, problem, algo) + @test_skip ncolors(result) >= ncolors(coloring(A, problem, optalgo)) + end +end + +@testset "Row coloring" begin + problem = ColoringProblem(; structure=:nonsymmetric, partition=:row) + for (m, n, p) in asymmetric_params + A = sprand(rng, m, n, p) + result = coloring(A, problem, algo) + @test_skip ncolors(result) >= ncolors(coloring(A, problem, optalgo)) + end +end + +@testset "Too big" begin + A = sprand(rng, Bool, 100, 100, 0.1) + optalgo_timelimit = OptimalColoringAlgorithm( + optimizer_with_attributes(HiGHS.Optimizer, "time_limit" => 10.0); # 1 second + silent=false, + assert_solved=false, + ) + @test_throws AssertionError coloring(A, ColoringProblem(), optalgo_timelimit) +end diff --git a/test/order.jl b/test/order.jl index 9bb9e6dc..9cae16f1 100644 --- a/test/order.jl +++ b/test/order.jl @@ -61,16 +61,14 @@ end; end; @testset "LargestFirst" begin - for has_diagonal in (false, true) - A = sparse([ - 0 1 0 0 - 1 0 1 1 - 0 1 0 1 - 0 1 1 0 - ]) - ag = AdjacencyGraph(A; has_diagonal) - @test vertices(ag, LargestFirst()) == [2, 3, 4, 1] - end + A = sparse([ + 0 1 0 0 + 1 0 1 1 + 0 1 0 1 + 0 1 1 0 + ]) + ag = AdjacencyGraph(A) + @test vertices(ag, LargestFirst()) == [2, 3, 4, 1] A = sparse([ 1 1 0 0 @@ -146,3 +144,99 @@ end; @test isperm(π) end end + +@testset "Multiple orders" begin + # I used brute force to find examples where LargestFirst is *strictly* better than NaturalOrder, just to check that the best order is indeed selected when multiple orders are provided + @testset "Column coloring" begin + A = [ + 0 0 1 1 + 0 1 0 1 + 0 0 1 1 + 1 1 0 0 + ] + problem = ColoringProblem{:nonsymmetric,:column}() + algo = GreedyColoringAlgorithm(NaturalOrder()) + better_algo = GreedyColoringAlgorithm((NaturalOrder(), LargestFirst())) + @test ncolors(coloring(A, problem, better_algo)) < + ncolors(coloring(A, problem, algo)) + end + @testset "Row coloring" begin + A = [ + 1 0 0 0 + 0 0 1 0 + 0 1 1 1 + 1 0 0 1 + ] + problem = ColoringProblem{:nonsymmetric,:row}() + algo = GreedyColoringAlgorithm(NaturalOrder()) + better_algo = GreedyColoringAlgorithm((NaturalOrder(), LargestFirst())) + @test ncolors(coloring(A, problem, better_algo)) < + ncolors(coloring(A, problem, algo)) + end + @testset "Star coloring" begin + A = [ + 0 1 0 1 1 + 1 1 0 1 0 + 0 0 1 0 1 + 1 1 0 1 0 + 1 0 1 0 0 + ] + problem = ColoringProblem{:symmetric,:column}() + algo = GreedyColoringAlgorithm(NaturalOrder()) + better_algo = GreedyColoringAlgorithm((NaturalOrder(), LargestFirst())) + @test ncolors(coloring(A, problem, better_algo)) < + ncolors(coloring(A, problem, algo)) + end + @testset "Acyclic coloring" begin + A = [ + 1 0 0 0 0 1 0 + 0 0 0 1 0 0 0 + 0 0 0 1 0 0 0 + 0 1 1 1 0 1 1 + 0 0 0 0 0 0 1 + 1 0 0 1 0 0 1 + 0 0 0 1 1 1 1 + ] + problem = ColoringProblem{:symmetric,:column}() + algo = GreedyColoringAlgorithm{:substitution}(NaturalOrder()) + better_algo = GreedyColoringAlgorithm{:substitution}(( + NaturalOrder(), LargestFirst() + )) + @test ncolors(coloring(A, problem, better_algo)) < + ncolors(coloring(A, problem, algo)) + end + @testset "Star bicoloring" begin + A = [ + 0 1 0 0 0 + 1 0 1 0 0 + 0 1 0 0 1 + 0 0 0 0 0 + 0 0 1 0 1 + ] + problem = ColoringProblem{:nonsymmetric,:bidirectional}() + algo = GreedyColoringAlgorithm(NaturalOrder()) + better_algo = GreedyColoringAlgorithm((NaturalOrder(), LargestFirst())) + @test ncolors(coloring(A, problem, better_algo)) < + ncolors(coloring(A, problem, algo)) + end + @testset "Acyclic bicoloring" begin + A = [ + 0 1 0 1 1 0 1 0 1 + 1 0 0 0 0 0 0 0 1 + 0 0 0 0 0 0 0 0 0 + 1 0 0 1 1 0 1 0 0 + 1 0 0 1 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 + 1 0 0 1 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 + 1 1 0 0 0 0 0 0 0 + ] + problem = ColoringProblem{:nonsymmetric,:bidirectional}() + algo = GreedyColoringAlgorithm{:substitution}(NaturalOrder()) + better_algo = GreedyColoringAlgorithm{:substitution}(( + NaturalOrder(), LargestFirst() + )) + @test ncolors(coloring(A, problem, better_algo)) < + ncolors(coloring(A, problem, algo)) + end +end diff --git a/test/random.jl b/test/random.jl index 02c1a3a1..406dfea9 100644 --- a/test/random.jl +++ b/test/random.jl @@ -81,10 +81,10 @@ end; problem = ColoringProblem(; structure=:nonsymmetric, partition=:bidirectional) @testset for algo in ( GreedyColoringAlgorithm( - RandomOrder(rng); postprocessing=false, decompression=:direct + RandomOrder(StableRNG(0), 0); postprocessing=false, decompression=:direct ), GreedyColoringAlgorithm( - RandomOrder(rng); postprocessing=true, decompression=:direct + RandomOrder(StableRNG(0), 0); postprocessing=true, decompression=:direct ), ) @testset "$((; m, n, p))" for (m, n, p) in asymmetric_params @@ -102,10 +102,10 @@ end; problem = ColoringProblem(; structure=:nonsymmetric, partition=:bidirectional) @testset for algo in ( GreedyColoringAlgorithm( - RandomOrder(rng); postprocessing=false, decompression=:substitution + RandomOrder(StableRNG(0), 0); postprocessing=false, decompression=:substitution ), GreedyColoringAlgorithm( - RandomOrder(rng); postprocessing=true, decompression=:substitution + RandomOrder(StableRNG(0), 0); postprocessing=true, decompression=:substitution ), ) @testset "$((; m, n, p))" for (m, n, p) in asymmetric_params diff --git a/test/result.jl b/test/result.jl index 0bea3353..49344aab 100644 --- a/test/result.jl +++ b/test/result.jl @@ -1,3 +1,4 @@ +using SparseMatrixColorings using SparseMatrixColorings: group_by_color, UnsupportedDecompressionError using Test @@ -20,7 +21,7 @@ using Test end @testset "Empty compression" begin - A = rand(10, 10) + A = zeros(Bool, 10, 10) color = zeros(Int, 10) problem = ColoringProblem{:nonsymmetric,:column}() algo = ConstantColoringAlgorithm(A, color; partition=:column) diff --git a/test/runtests.jl b/test/runtests.jl index 81c4ecc9..ec679d4e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,7 +1,6 @@ using Aqua using Documenter using JET -using JuliaFormatter using SparseMatrixColorings using Test @@ -13,20 +12,17 @@ include("utils.jl") @testset verbose = true "SparseMatrixColorings" begin if get(ENV, "JULIA_SMC_TEST_GROUP", nothing) == "GPU" @testset "CUDA" begin - using CUDA + using CUDA, cuSPARSE include("cuda.jl") end else @testset verbose = true "Code quality" begin @testset "Aqua" begin - Aqua.test_all(SparseMatrixColorings; stale_deps=(; ignore=[:Requires],)) + Aqua.test_all(SparseMatrixColorings; undocumented_names=true) end @testset "JET" begin - JET.test_package(SparseMatrixColorings; target_defined_modules=true) - end - @testset "JuliaFormatter" begin - @test JuliaFormatter.format( - SparseMatrixColorings; verbose=false, overwrite=false + JET.test_package( + SparseMatrixColorings; target_modules=(SparseMatrixColorings,) ) end @testset "Doctests" begin @@ -58,6 +54,9 @@ include("utils.jl") @testset "Constant coloring" begin include("constant.jl") end + @testset "Optimal coloring" begin + include("optimal.jl") + end @testset "ADTypes coloring algorithms" begin include("adtypes.jl") end @@ -84,7 +83,10 @@ include("utils.jl") end @testset verbose = true "Performance" begin @testset "Type stability" begin - include("type_stability.jl") + if VERSION < v"1.12" + # TODO: fix JET misbehaving + include("type_stability.jl") + end end @testset "Allocations" begin include("allocations.jl") diff --git a/test/structured.jl b/test/structured.jl index 23ba967f..857ebbfc 100644 --- a/test/structured.jl +++ b/test/structured.jl @@ -2,57 +2,86 @@ using ArrayInterface: ArrayInterface using BandedMatrices: BandedMatrix, brand using BlockBandedMatrices: BandedBlockBandedMatrix, BlockBandedMatrix using LinearAlgebra +using SparseArrays using SparseMatrixColorings using Test @testset "Diagonal" begin - for n in (1, 2, 10, 100) - A = Diagonal(rand(n)) - test_structured_coloring_decompression(A) + @testset for algo in [GreedyColoringAlgorithm(), StructuredColoringAlgorithm()] + for n in (1, 2, 10, 100) + A = Diagonal(rand(n)) + test_structured_coloring_decompression(A, algo) + end end end; @testset "Bidiagonal" begin - for n in (2, 10, 100) - A1 = Bidiagonal(rand(n), rand(n - 1), :U) - A2 = Bidiagonal(rand(n), rand(n - 1), :L) - test_structured_coloring_decompression(A1) - test_structured_coloring_decompression(A2) + @testset for algo in [GreedyColoringAlgorithm(), StructuredColoringAlgorithm()] + for n in (2, 10, 100) + A1 = Bidiagonal(rand(n), rand(n - 1), :U) + A2 = Bidiagonal(rand(n), rand(n - 1), :L) + test_structured_coloring_decompression(A1, algo) + test_structured_coloring_decompression(A2, algo) + end end end; @testset "Tridiagonal" begin - for n in (2, 10, 100) - A = Tridiagonal(rand(n - 1), rand(n), rand(n - 1)) - test_structured_coloring_decompression(A) + @testset for algo in [GreedyColoringAlgorithm(), StructuredColoringAlgorithm()] + for n in (2, 10, 100) + A = Tridiagonal(rand(n - 1), rand(n), rand(n - 1)) + test_structured_coloring_decompression(A, algo) + end end end; @testset "BandedMatrices" begin - @testset for (m, n) in [(10, 20), (20, 10)], l in 0:5, u in 0:5 - A = brand(m, n, l, u) - test_structured_coloring_decompression(A) + @testset for algo in [GreedyColoringAlgorithm(), StructuredColoringAlgorithm()] + @testset for (m, n) in [(10, 20), (20, 10)], l in 0:5, u in 0:5 + A = brand(m, n, l, u) + test_structured_coloring_decompression(A, algo) + end end end; @testset "BlockBandedMatrices" begin - for (mb, nb) in [(10, 20), (20, 10)], lb in 0:3, ub in 0:3, _ in 1:10 - rows = rand(1:5, mb) - cols = rand(1:5, nb) - A = BlockBandedMatrix{Float64}(rand(sum(rows), sum(cols)), rows, cols, (lb, ub)) - test_structured_coloring_decompression(A) + @testset for algo in [GreedyColoringAlgorithm(), StructuredColoringAlgorithm()] + for (mb, nb) in [(10, 20), (20, 10)], lb in 0:3, ub in 0:3, _ in 1:10 + rows = rand(1:5, mb) + cols = rand(1:5, nb) + A = BlockBandedMatrix{Float64}(rand(sum(rows), sum(cols)), rows, cols, (lb, ub)) + test_structured_coloring_decompression(A, algo) + end end end; @testset "BandedBlockBandedMatrices" begin - for (mb, nb) in [(10, 20), (20, 10)], lb in 0:3, ub in 0:3, _ in 1:10 - rows = rand(5:10, mb) - cols = rand(5:10, nb) - λ = rand(0:5) - μ = rand(0:5) - A = BandedBlockBandedMatrix{Float64}( - rand(sum(rows), sum(cols)), rows, cols, (lb, ub), (λ, μ) - ) - test_structured_coloring_decompression(A) + @testset for algo in [GreedyColoringAlgorithm(), StructuredColoringAlgorithm()] + for (mb, nb) in [(10, 20), (20, 10)], lb in 0:3, ub in 0:3, _ in 1:10 + rows = rand(5:10, mb) + cols = rand(5:10, nb) + λ = rand(0:5) + μ = rand(0:5) + A = BandedBlockBandedMatrix{Float64}( + rand(sum(rows), sum(cols)), rows, cols, (lb, ub), (λ, μ) + ) + test_structured_coloring_decompression(A, algo) + end end end; + +# See https://github.com/JuliaDiff/SparseMatrixColorings.jl/pull/299 +@testset "SparsityPatternCSC $T" for T in [Int, Float32] + S = sparse(T[ + 0 0 1 1 0 1 + 1 0 0 0 1 0 + 0 1 0 0 1 0 + 0 1 1 0 0 0 + ]) + P = SparseMatrixColorings.SparsityPatternCSC(S) + problem = ColoringProblem() + algo = GreedyColoringAlgorithm() + result = coloring(P, problem, algo) + B = compress(S, result) + @test decompress(B, result) isa SparseMatrixCSC{T,Int} +end; diff --git a/test/type_stability.jl b/test/type_stability.jl index 5a488596..3c6f9f02 100644 --- a/test/type_stability.jl +++ b/test/type_stability.jl @@ -40,11 +40,21 @@ rng = StableRNG(63) ColoringProblem(; structure, partition), GreedyColoringAlgorithm(order; decompression), ) + @test_opt coloring( + A, + ColoringProblem(; structure, partition), + GreedyColoringAlgorithm((NaturalOrder(), order); decompression), + ) @inferred coloring( A, ColoringProblem(; structure, partition), GreedyColoringAlgorithm(order; decompression), ) + @inferred coloring( + A, + ColoringProblem(; structure, partition), + GreedyColoringAlgorithm((NaturalOrder(), order); decompression), + ) end end end; diff --git a/test/utils.jl b/test/utils.jl index bb80f95f..15aa3f9f 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -12,9 +12,16 @@ using SparseMatrixColorings: respectful_similar, structurally_orthogonal_columns, symmetrically_orthogonal_columns, - structurally_biorthogonal + structurally_biorthogonal, + substitutable_columns, + substitutable_bidirectional, + rank_nonzeros_from_trees using Test +const _ALL_ORDERS = ( + NaturalOrder(), LargestFirst(), SmallestLast(), IncidenceDegree(), DynamicLargestFirst() +) + function test_coloring_decompression( A0::AbstractMatrix, problem::ColoringProblem{structure,partition}, @@ -73,7 +80,6 @@ function test_coloring_decompression( end @testset "Recoverability" begin - # TODO: find tests for recoverability for substitution decompression if decompression == :direct if structure == :nonsymmetric if partition == :column @@ -93,6 +99,15 @@ function test_coloring_decompression( end end + @testset "Substitutable" begin + if decompression == :substitution + if structure == :symmetric + rank_nonzeros = rank_nonzeros_from_trees(result) + @test substitutable_columns(A0, rank_nonzeros, color) + end + end + end + @testset "Single-color decompression" begin if decompression == :direct # TODO: implement for :substitution too A2 = respectful_similar(A, eltype(B)) @@ -169,6 +184,22 @@ function test_coloring_decompression( @show color_vec end end + + @testset "More orders is better" begin + more_orders = (algo.orders..., _ALL_ORDERS...) + better_algo = GreedyColoringAlgorithm{decompression}( + more_orders; algo.postprocessing + ) + all_algos = [ + GreedyColoringAlgorithm{decompression}(order; algo.postprocessing) for + order in more_orders + ] + result = coloring(A0, problem, algo) + better_result = coloring(A0, problem, better_algo) + all_results = [coloring(A0, problem, _algo) for _algo in all_algos] + @test ncolors(better_result) <= ncolors(result) + @test ncolors(better_result) == minimum(ncolors, all_results) + end end function test_bicoloring_decompression( @@ -215,13 +246,39 @@ function test_bicoloring_decompression( @test structurally_biorthogonal(A0, row_color, column_color) end end + + if decompression == :substitution + @testset "Substitutable" begin + rank_nonzeros = rank_nonzeros_from_trees(result) + @test substitutable_bidirectional( + A0, rank_nonzeros, row_color, column_color + ) + end + end + end + + @testset "More orders is better" begin + more_orders = (algo.orders..., _ALL_ORDERS...) + better_algo = GreedyColoringAlgorithm{decompression}( + more_orders; algo.postprocessing + ) + all_algos = [ + GreedyColoringAlgorithm{decompression}(order; algo.postprocessing) for + order in more_orders + ] + result = coloring(A0, problem, algo) + better_result = coloring(A0, problem, better_algo) + all_results = [coloring(A0, problem, _algo) for _algo in all_algos] + @test ncolors(better_result) <= ncolors(result) + @test ncolors(better_result) == minimum(ncolors, all_results) end end -function test_structured_coloring_decompression(A::AbstractMatrix) +function test_structured_coloring_decompression( + A::AbstractMatrix, algo=StructuredColoringAlgorithm() +) column_problem = ColoringProblem(; structure=:nonsymmetric, partition=:column) row_problem = ColoringProblem(; structure=:nonsymmetric, partition=:row) - algo = GreedyColoringAlgorithm() # Column result = coloring(A, column_problem, algo) @@ -231,9 +288,8 @@ function test_structured_coloring_decompression(A::AbstractMatrix) @test D == A @test nameof(typeof(D)) == nameof(typeof(A)) @test structurally_orthogonal_columns(A, color) - if VERSION >= v"1.10" || A isa Union{Diagonal,Bidiagonal,Tridiagonal} - # banded matrices not supported by ArrayInterface on Julia 1.6 - # @test color == ArrayInterface.matrix_colors(A) # TODO: uncomment + if algo isa StructuredColoringAlgorithm + @test color == ArrayInterface.matrix_colors(A) end # Row