julialinearalgebra / hierarchicalmatrices.jl

Julia package for hierarchical matrices

HierarchicalMatrices.jl

This package provides a flexible framework for hierarchical data types in Julia.

Create your own hierarchical matrix as simply as:

julia> using HierarchicalMatrices

julia> @hierarchical MyHierarchicalMatrix LowRankMatrix Matrix

The invocation of the @hierarchical macro creates an abstract supertype AbstractMyHierarchicalMatrix{T} <: AbstractMatrix{T} and the immutable type MyHierarchicalMatrix, endowing it with fields of HierarchicalMatrixblocks, LowRankMatrixblocks, Matrixblocks, and a matrix of integers, assigned, to determine which type of block is active. The package comes pre-loaded with a HierarchicalMatrix.

See the example on speeding up the matrix-vector product with Cauchy matrices.

Implementation

A straightforward implementation of hierarchical (self-referential) data types might suffer from Russell's paradox. In the context of types, Russell's paradox states that either you know the type, or you know its fields, but neither concretely. On one side of the paradox, you end up with type-stable constructors and type-unstable getters; on the other side, you are stuck with type-unstable constructors and type-stable getters.

This implementation of hierarchical data types avoids Russell's paradox at the cost of restricting the entire list of concrete matrix types that are the fields upon construction. This allows for fast and type-stable setters and getters. Enjoy!

hierarchicalmatrices.jl's Issues

Error tagging new release

The tag name "0.1.1" is not of the appropriate SemVer form (vX.Y.Z).
cc: @MikaelSlevinsky

Error tagging new release

The URL of this package does not match that stored in METADATA.jl.
cc: @MikaelSlevinsky

Error tagging new release

The URL of this package does not match that stored in METADATA.jl.
cc: @MikaelSlevinsky

Use femtocleaner bot to update deprecated syntax

https://github.com/apps/femtocleaner

HierarchicalMatrices should use BlockArrays

BlockArrays has a well-developed block matrix structure, so this would avoid needing your own Block.

KernelMatrix goes into infinite loop

using LinearAlgebra
using HierarchicalMatrices

h = 0.05
L = 5.0
N = Int(floor(L/h))
xs = range(-L,L,length=N)

K(x,y) = 1/(1+abs(x-y))


KM = KernelMatrix(K,Array(xs),Array(xs),-L,L,-L,L)

goes into an infinite loop. h=0.1 works.

Move LowRankMatrix to LowRankApprox (or new package LowRankMatrices.jl)

LowRankMatrix is not actually a HierarchicalMatrix, so maybe it lives better in LowRankApprox. This would also make it easier to use the LowRankApprox algorithms.

How best to use shared memory parallelism in Julia?

A ccall acts as a thread barrier, so even if it's calling BLAS which is multithreaded, the overhead might be high for exceptionally low-rank data structures. Perhaps there's a better way to thread from the top down.

Different name for partial-update functions?

The signatures of A_mul_B! and scale! and similar here that take start indices that don't match the signatures used by those functions in base don't really need to use the same name here, do they? It might be clearer to distinguish them by more than just the number of inputs, since Base is never going to lower to or call these signatures.

norm(matrix) -> opnorm(matrix) in Julia 0.7

It looks like may be using norm(matrix). In Julia 0.7, this will compute the Frobenius norm (vecnorm in Julia 0.6), due to JuliaLang/julia#27401. If you want the induced/operator norm as in Julia 0.6, use opnorm(matrix) instead, or Compat.opnorm(matrix) to work in 0.6 and 0.7 (JuliaLang/Compat.jl#577).

Note that, for testing purposes, rather than @test norm(A - B) ≤ tol, it is usually preferred to do @test A ≈ B or @test A ≈ B rtol=... (which uses isapprox).

Error tagging new release

The URL of this package does not match that stored in METADATA.jl.
cc: @MikaelSlevinsky

Move HierarchicalVector and hierarchicalsolve here?

Is there a reason these are still in SingularIntegralEquations.jl?

2D with reshaped array

Really nice work! I am wondering if this is also suitable for the convolution of a 2D field? Do you require the indices of x be related to the physical position it represents? Thanks!

Ke

Move to JuliaMatrices?

The benefit is that members (e.g. me and @andreasnoack) can help maintain the package, incorporate apps like Femtobots, etc.