Repository navigation
feature/harithmetic - #36
Merged
Merged
Conversation
…s keyword
HMatrix no longer stores its near and far field as BlockSparseMatrices
objects. The new HMatrixCore (HMatrix remains as an alias) keeps the block
partition itself, which the hierarchical arithmetic builds on:
- the retained block tree, the space ordering and the assembly tolerance `tol`
- dense near blocks (Vector{Matrix{T}}) and low-rank far blocks
(Vector{LowRankMatrix{T}}), each stored once
- a canonical CSR ordering of the terminal blocks and a CSC traversal map
that refers back to the CSR entries without duplicating blocks
Assembly
- The tree backend reports the terminal block partition
(terminalinteractions, implemented in ACAH2Trees); near blocks are filled
through assemble_blocks and far blocks are compressed in parallel.
- Transpose-symmetric assembly stores each unordered pair once and refers to
it transposed for the opposite block. It now also supports
PermuteSpaceInPlace.
- permutation(tree) only reads the tree. PermuteSpaceInPlace renumbers the
tree with the new _settreeorder! after permuting the space, so querying the
permutation (e.g. when building a workspace) no longer alters the tree.
- LowRankMatrix is mutable and no longer carries a product buffer.
- copy(H) copies the numerical blocks and shares tree and block index.
- nearmatrix/farmatrix return the near or far part of the stored matrix;
nnz/storage count the new storage.
Matrix-vector products
- MatVecWorkspace(H) holds the matrix, a precomputed schedule of
cost-balanced output chunks and all buffers, and is used in place of H:
products, transpose and adjoint products, mul! with α and β, several
right-hand sides, real matrices with complex vectors. Repeated products do
not allocate numerical storage.
- Products with H itself build a temporary workspace on every call and warn
once.
Threading interface (breaking)
- The `scheduler`, `mvpscheduler` and `nchunks` keywords are replaced by
`ntasks` (default Threads.nthreads()) in HMatrix, assemble, assemble_blocks,
assemblenears, assemblefars and MatVecWorkspace. The package picks the
scheduling: a greedy task pool for assembly, eight cost-balanced chunks per
task for products, serial execution for ntasks=1.
- assemble_blocks, assemblenears and assemblefars now default to
Threads.nthreads() tasks instead of serial execution. Callers passing
`scheduler=` (NestedCrossApproximation) have to switch to `ntasks`.
BlockSparseMatrices is still used by the standalone assemblenears and
assemblefars.
Tests cover the new storage, the tree index and the products; the symmetric
assembly tests are rewritten for the new storage.
H-matrix arithmetic through Julia's standard functions, for matrices with full (non-symmetric) storage: - C = A + B, A - B, α*A, axpy!, axpby!, rmul!, lmul! - C = A * B and mul!(C, A, B, α, β); sums need identical block partitions, products matching cluster trees - recompress!(H): truncates every far block to H.tol - lu(H), lu!(H) and solves with vectors and dense matrices: \, ldiv!, also with transpose and adjoint factors; issuccess The tolerance chosen at assembly governs every later recompression, relative to each block in the spectral norm; for operands with different tolerances the larger one is used. Workspaces Each operation has an internal, reusable workspace that holds its plan and buffers and is passed in place of the matrix: AdditionWorkspace, ProductWorkspace, RecompressionWorkspace, LUWorkspace and SolveWorkspace. LUWorkspace is bound to the partition, not to one matrix, so a frequency sweep on one geometry factorizes every matrix without planning again. Operations called without a workspace build a temporary one. Products - Low-rank terminals multiply whole hierarchical panels; a low-rank terminal above several destinations is formed and recompressed once on a virtual aggregation node. - Cheap far destinations are recompressed directly from the implicit operator C + α*A*B; remaining updates are batched or summed into a dense core when that is cheaper. - Wide recompressions use a randomized range finder with an a posteriori error bound followed by a small SVD. The sketches use a fixed-seed RNG (new dependency on Random), so results are reproducible. LU factorization - Left-looking hierarchical LU, P*A ≈ L*U, with pivoting restricted to each dense diagonal leaf. A zero pivot throws LeafLUFailure, or gives issuccess(F) == false with check=false. The factorization (HMatrixLU) reports the local recompression status (converged, rank_limited, max_block_error, flushes). Parallelism - LU and products build a data-flow task graph once in their workspace: a task depends on the last writer of every block it reads or writes and on the earlier readers of the blocks it writes. ntasks workers then take the earliest ready task. The results are bitwise identical for every ntasks. - Sums run over independent blocks; recompression on ntasks tasks with one scratch per task. - Measured with one BLAS thread: 30k unknowns 5.3 s serial, 0.55 s on 16 tasks; 100k unknowns 29 s serial, 2.7 s on 32 tasks. Beyond about 16 tasks, process-level overhead (OpenBLAS buffer lock, kernel time) limits the scaling on the test machine. example/pointcloudlu_benchmark.jl compares the LU with HMatrices.jl on identical partitions (50k unknowns: 10.8 s vs 12.2 s on a box partition). Tests for every operation, including parallel-versus-serial equality of LU, products and sums. Twenty assertions in test_hierarchicalmultiplication.jl still expect the earlier buffer sizes, pointers and flush counts and fail; they are left for the upcoming cleanup.
- Terminal blocks of a block tree in compressed row storage (interactions, symmetricinteractions); ACAH2Trees provides the cluster interface and the traversals. Rename _tree to buildtree. - Separate near- and far-field assembly into one HMatrix with BlockCSR/BlockCSC block storage; drop the BlockSparseMatrices dependency. - Move the hierarchical arithmetic and the matrix-vector product to src/hmatrix/arithmetic, each file starting with its workspace. Split the LU solves into solve.jl and add H \ b. - Rewrite the parallel matrix-vector product; transposed products avoid the OpenBLAS scratch-buffer lock. - Reorganize the test suite by source directory.
- Examples assemble the operator and solve with \; add the point-cloud example and the BEAST adapter for combined operators. - Document the arithmetic, LU and solves in the manual and the details section.
…ace, CompScienceMeshes 0.13
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
No description provided.