From 1aec62a3d6e21c52959f900a3124ca87b43b4968 Mon Sep 17 00:00:00 2001 From: aliat Date: Wed, 13 Aug 2025 17:36:14 +0200 Subject: [PATCH 1/4] different pivoting for hmatrices --- ext/ACAHMatrices/hmatrices_interface.jl | 87 ++++++++++++++++++++++++- src/convergence/randomsampling.jl | 4 +- test/test_ACAHMatrices.jl | 56 +++++++++++++++- 3 files changed, 138 insertions(+), 9 deletions(-) diff --git a/ext/ACAHMatrices/hmatrices_interface.jl b/ext/ACAHMatrices/hmatrices_interface.jl index 33bc612..49d0bc2 100644 --- a/ext/ACAHMatrices/hmatrices_interface.jl +++ b/ext/ACAHMatrices/hmatrices_interface.jl @@ -18,12 +18,13 @@ function (aca::ACA)( colBuffer = reshape( bufs.B.data[1:(length(irange) * maxrank)], (length(irange), maxrank) ) - + update_conv_idcs!(aca.convergence, K, irange, jrange) cur_aca = AdaptiveCrossApproximation.ACA(; - rowpivoting=AdaptiveCrossApproximation.MaximumValue(zeros(Bool, length(irange))), - columnpivoting=AdaptiveCrossApproximation.MaximumValue(zeros(Bool, length(jrange))), + rowpivoting=update_pivot_idcs(K, aca.rowpivoting, irange, jrange, true), + columnpivoting=update_pivot_idcs(K, aca.columnpivoting, irange, jrange, false), convergence=aca.convergence, ) + npivots, U, V = cur_aca( K, rowBuffer, colBuffer, maxrank; rowidcs=irange, colidcs=jrange ) @@ -32,6 +33,56 @@ function (aca::ACA)( ) end +function update_pivot_idcs( + _::HMatrices.PermutedMatrix{TT,T}, + pivoting::AdaptiveCrossApproximation.ValuePivStrat, + rowidcs::UnitRange{Int}, + colidcs::UnitRange{Int}, + is_row::Bool, +) where {TT<:HMatrices.KernelMatrix,T} + len = is_row ? length(rowidcs) : length(colidcs) + return pivoting(len) +end + +function update_pivot_idcs( + _::HMatrices.PermutedMatrix{TT,T}, + pivoting::AdaptiveCrossApproximation.GeoPivStrat, + rowidcs::UnitRange{Int}, + colidcs::UnitRange{Int}, + is_row::Bool, +) where {TT<:HMatrices.KernelMatrix,T} + idcsvec = is_row ? Vector(rowidcs) : Vector(colidcs) + return pivoting(idcsvec) +end + +function update_pivot_idcs( + K::HMatrices.PermutedMatrix{TT,T}, + pivoting::AdaptiveCrossApproximation.ConvPivStrat, + rowidcs::UnitRange{Int}, + colidcs::UnitRange{Int}, + _::Bool, +) where {TT<:HMatrices.KernelMatrix,T} + update_conv_idcs!(pivoting.convcrit, K, rowidcs, colidcs) + return typeof(pivoting)(pivoting.convcrit, pivoting.rc) +end + +function update_pivot_idcs( + K::HMatrices.PermutedMatrix{TT,T}, + pivoting::AdaptiveCrossApproximation.CombinedPivStrat, + rowidcs::UnitRange{Int}, + colidcs::UnitRange{Int}, + is_row::Bool, +) where {TT<:HMatrices.KernelMatrix,T} + curr_strats = Vector{AdaptiveCrossApproximation.PivStrat}( + undef, length(pivoting.strats) + ) + for (i, strat) in enumerate(pivoting.strats) + curr_strats[i] = update_pivot_idcs(K, strat, rowidcs, colidcs, is_row) + end + update_conv_idcs!(pivoting.convcrit, K, rowidcs, colidcs) + return AdaptiveCrossApproximation.CombinedPivStrat(pivoting.convcrit, curr_strats) +end + function Base.resize!(A::HMatrices.VectorOfVectors, m::Int, n::Int) A.m = m ie = m * n @@ -41,6 +92,36 @@ function Base.resize!(A::HMatrices.VectorOfVectors, m::Int, n::Int) return A.k = n end +function update_conv_idcs!( + convcrit::AdaptiveCrossApproximation.CombinedConvCrit, + K::HMatrices.PermutedMatrix{TT,T}, + irange::UnitRange{Int}, + jrange::UnitRange{Int}, +) where {TT<:HMatrices.KernelMatrix,T} + for crit in convcrit.crits + update_conv_idcs!(crit, K, irange, jrange) + end +end + +update_conv_idcs!( + convcrit::AdaptiveCrossApproximation.ConvCrit, + K::HMatrices.PermutedMatrix{TT,T}, + irange::UnitRange{Int}, + jrange::UnitRange{Int}, +) where {TT<:HMatrices.KernelMatrix,T} = nothing + +function update_conv_idcs!( + convcrit::AdaptiveCrossApproximation.RandomSampling{F,G}, + K::HMatrices.PermutedMatrix{TT,T}, + irange::UnitRange{Int}, + jrange::UnitRange{Int}, +) where {F<:Real,G,TT<:HMatrices.KernelMatrix,T} + convcrit.indices = hcat( + rand(1:length(irange), convcrit.nsamples), rand(1:length(jrange), convcrit.nsamples) + ) + return convcrit.rest = [K.data[rc[1], rc[2]][1] for rc in eachrow(convcrit.indices)] +end + Base.size(K::HMatrices.VectorOfVectors) = K.m, K.k function AdaptiveCrossApproximation.nextrc!( diff --git a/src/convergence/randomsampling.jl b/src/convergence/randomsampling.jl index d543a72..f7830a2 100644 --- a/src/convergence/randomsampling.jl +++ b/src/convergence/randomsampling.jl @@ -12,9 +12,7 @@ end function RandomSampling( ::Type{K}; factor::F=1.0, nsamples::Int=0, tol::F=1e-4 ) where {K,F<:Real} - return RandomSampling{F,K}( - F(0.0), nsamples, factor, zeros(Int, 0, 0), zeros(K, 0, 0), tol - ) + return RandomSampling{F,K}(F(0.0), nsamples, factor, zeros(Int, 0, 0), zeros(K, 0), tol) end tolerance(cc::RandomSampling) = cc.tol diff --git a/test/test_ACAHMatrices.jl b/test/test_ACAHMatrices.jl index 844df8e..c01bb8c 100644 --- a/test/test_ACAHMatrices.jl +++ b/test/test_ACAHMatrices.jl @@ -2,6 +2,7 @@ using BEAST using HMatrices using CompScienceMeshes using AdaptiveCrossApproximation +using Random using Test using LinearAlgebra @@ -13,18 +14,67 @@ xclt = ClusterTree(spaceX.pos) K = HMatrices.KernelMatrix(op, spaceX, spaceX); -rtols = [10.0^i for i in collect(-4:-1:-10)] +rtols = [1e-2, 1e-4, 1e-6, 1e-8, 1e-10, 1e-12, 1e-14] tst_vec = rand(dim) fullmat = BEAST.assemble(op, spaceX, spaceX) trueResult = fullmat * tst_vec; -## +# PivStrat :: MaximumValue for rtol in rtols - K = HMatrices.KernelMatrix(op, spaceX, spaceX) aca = ACA(; convergence=FNormEstimator(0.0, rtol)) hmat = HMatrices.assemble_hmatrix(K; comp=aca) @test size(hmat, 1) == size(fullmat, 1) @test size(hmat, 2) == size(fullmat, 2) @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol end +# PivStrat :: leja2 +for rtol in rtols + aca = ACA(; rowpivoting=Leja2(spaceX.pos), convergence=FNormEstimator(0.0, rtol)) + hmat = HMatrices.assemble_hmatrix(K; comp=aca) + @test size(hmat, 1) == size(fullmat, 1) + @test size(hmat, 2) == size(fullmat, 2) + @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol + + aca = ACA(; columnpivoting=Leja2(spaceX.pos), convergence=FNormEstimator(0.0, rtol)) + hmat = HMatrices.assemble_hmatrix(K; comp=aca) + @test size(hmat, 1) == size(fullmat, 1) + @test size(hmat, 2) == size(fullmat, 2) + @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol +end + +# PivStrat :: FillDistance +for rtol in rtols + aca = ACA(; rowpivoting=FillDistance(spaceX.pos), convergence=FNormEstimator(0.0, rtol)) + hmat = HMatrices.assemble_hmatrix(K; comp=aca) + @test size(hmat, 1) == size(fullmat, 1) + @test size(hmat, 2) == size(fullmat, 2) + @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol + + aca = ACA(; + columnpivoting=FillDistance(spaceX.pos), convergence=FNormEstimator(0.0, rtol) + ) + hmat = HMatrices.assemble_hmatrix(K; comp=aca) + @test size(hmat, 1) == size(fullmat, 1) + @test size(hmat, 2) == size(fullmat, 2) + @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol +end + +# PivStrat :: CombinedPivStrat +for rtol in [1e-2, 1e-4, 1e-6, 1e-8, 1e-10]#, 1e-12, 1e-14] + Random.seed!(1) + cc1 = AdaptiveCrossApproximation.FNormEstimator(rtol) + cc2 = AdaptiveCrossApproximation.RandomSampling( + Float64; nsamples=100, factor=1.0, tol=rtol + ) + convergence = AdaptiveCrossApproximation.CombinedConvCrit([cc1, cc2], zeros(Bool, 2)) + + ps1 = MaximumValue() + ps2 = AdaptiveCrossApproximation.RandomSamplingPivoting(cc2, 1) + rp = AdaptiveCrossApproximation.CombinedPivStrat(convergence, [ps1, ps2]) + aca = ACA(; rowpivoting=rp, convergence=convergence) + hmat = HMatrices.assemble_hmatrix(K; comp=aca) + @test size(hmat, 1) == size(fullmat, 1) + @test size(hmat, 2) == size(fullmat, 2) + @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol +end From 78abe253043a6a6bad07f2e98cfbe82993b431f8 Mon Sep 17 00:00:00 2001 From: aliat Date: Thu, 14 Aug 2025 18:39:55 +0200 Subject: [PATCH 2/4] different pivot_linting fail --- Project.toml | 5 +- ext/ACAHMatrices/hmatrices_interface.jl | 91 +------------------- src/aca.jl | 60 ++++++++++++- src/convergence/combinedconvcrit.jl | 14 +++ src/convergence/estimation.jl | 5 +- src/convergence/extrapolation.jl | 4 +- src/convergence/randomsampling.jl | 6 ++ src/pivoting/combinedpivstrat.jl | 14 ++- src/pivoting/maxvalue.jl | 1 + src/pivoting/randomsampling.jl | 4 + test/runtests.jl | 3 +- test/test_ACAHMatrices.jl | 109 +++++++++++++----------- 12 files changed, 167 insertions(+), 149 deletions(-) diff --git a/Project.toml b/Project.toml index 4b0b28d..f8e3b78 100644 --- a/Project.toml +++ b/Project.toml @@ -4,13 +4,11 @@ authors = ["JoshuaTetzner "] version = "0.1.0" [deps] +HMatrices = "8646bddf-ab1c-4fa7-9c51-ba187d647618" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" -[weakdeps] -HMatrices = "8646bddf-ab1c-4fa7-9c51-ba187d647618" - [extensions] ACAHMatrices = ["HMatrices"] @@ -22,6 +20,7 @@ StaticArrays = "1.9.13" [extras] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +HMatrices = "8646bddf-ab1c-4fa7-9c51-ba187d647618" JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" diff --git a/ext/ACAHMatrices/hmatrices_interface.jl b/ext/ACAHMatrices/hmatrices_interface.jl index 49d0bc2..f10e659 100644 --- a/ext/ACAHMatrices/hmatrices_interface.jl +++ b/ext/ACAHMatrices/hmatrices_interface.jl @@ -18,71 +18,14 @@ function (aca::ACA)( colBuffer = reshape( bufs.B.data[1:(length(irange) * maxrank)], (length(irange), maxrank) ) - update_conv_idcs!(aca.convergence, K, irange, jrange) - cur_aca = AdaptiveCrossApproximation.ACA(; - rowpivoting=update_pivot_idcs(K, aca.rowpivoting, irange, jrange, true), - columnpivoting=update_pivot_idcs(K, aca.columnpivoting, irange, jrange, false), - convergence=aca.convergence, - ) - npivots, U, V = cur_aca( - K, rowBuffer, colBuffer, maxrank; rowidcs=irange, colidcs=jrange - ) + aca = aca(K, Vector(irange), Vector(jrange)) + npivots, U, V = aca(K, rowBuffer, colBuffer, maxrank; rowidcs=irange, colidcs=jrange) return HMatrices.RkMatrix( colBuffer[:, 1:npivots], Matrix(transpose(rowBuffer[1:npivots, :])) ) end -function update_pivot_idcs( - _::HMatrices.PermutedMatrix{TT,T}, - pivoting::AdaptiveCrossApproximation.ValuePivStrat, - rowidcs::UnitRange{Int}, - colidcs::UnitRange{Int}, - is_row::Bool, -) where {TT<:HMatrices.KernelMatrix,T} - len = is_row ? length(rowidcs) : length(colidcs) - return pivoting(len) -end - -function update_pivot_idcs( - _::HMatrices.PermutedMatrix{TT,T}, - pivoting::AdaptiveCrossApproximation.GeoPivStrat, - rowidcs::UnitRange{Int}, - colidcs::UnitRange{Int}, - is_row::Bool, -) where {TT<:HMatrices.KernelMatrix,T} - idcsvec = is_row ? Vector(rowidcs) : Vector(colidcs) - return pivoting(idcsvec) -end - -function update_pivot_idcs( - K::HMatrices.PermutedMatrix{TT,T}, - pivoting::AdaptiveCrossApproximation.ConvPivStrat, - rowidcs::UnitRange{Int}, - colidcs::UnitRange{Int}, - _::Bool, -) where {TT<:HMatrices.KernelMatrix,T} - update_conv_idcs!(pivoting.convcrit, K, rowidcs, colidcs) - return typeof(pivoting)(pivoting.convcrit, pivoting.rc) -end - -function update_pivot_idcs( - K::HMatrices.PermutedMatrix{TT,T}, - pivoting::AdaptiveCrossApproximation.CombinedPivStrat, - rowidcs::UnitRange{Int}, - colidcs::UnitRange{Int}, - is_row::Bool, -) where {TT<:HMatrices.KernelMatrix,T} - curr_strats = Vector{AdaptiveCrossApproximation.PivStrat}( - undef, length(pivoting.strats) - ) - for (i, strat) in enumerate(pivoting.strats) - curr_strats[i] = update_pivot_idcs(K, strat, rowidcs, colidcs, is_row) - end - update_conv_idcs!(pivoting.convcrit, K, rowidcs, colidcs) - return AdaptiveCrossApproximation.CombinedPivStrat(pivoting.convcrit, curr_strats) -end - function Base.resize!(A::HMatrices.VectorOfVectors, m::Int, n::Int) A.m = m ie = m * n @@ -92,36 +35,6 @@ function Base.resize!(A::HMatrices.VectorOfVectors, m::Int, n::Int) return A.k = n end -function update_conv_idcs!( - convcrit::AdaptiveCrossApproximation.CombinedConvCrit, - K::HMatrices.PermutedMatrix{TT,T}, - irange::UnitRange{Int}, - jrange::UnitRange{Int}, -) where {TT<:HMatrices.KernelMatrix,T} - for crit in convcrit.crits - update_conv_idcs!(crit, K, irange, jrange) - end -end - -update_conv_idcs!( - convcrit::AdaptiveCrossApproximation.ConvCrit, - K::HMatrices.PermutedMatrix{TT,T}, - irange::UnitRange{Int}, - jrange::UnitRange{Int}, -) where {TT<:HMatrices.KernelMatrix,T} = nothing - -function update_conv_idcs!( - convcrit::AdaptiveCrossApproximation.RandomSampling{F,G}, - K::HMatrices.PermutedMatrix{TT,T}, - irange::UnitRange{Int}, - jrange::UnitRange{Int}, -) where {F<:Real,G,TT<:HMatrices.KernelMatrix,T} - convcrit.indices = hcat( - rand(1:length(irange), convcrit.nsamples), rand(1:length(jrange), convcrit.nsamples) - ) - return convcrit.rest = [K.data[rc[1], rc[2]][1] for rc in eachrow(convcrit.indices)] -end - Base.size(K::HMatrices.VectorOfVectors) = K.m, K.k function AdaptiveCrossApproximation.nextrc!( diff --git a/src/aca.jl b/src/aca.jl index b801ea6..93f810f 100644 --- a/src/aca.jl +++ b/src/aca.jl @@ -20,13 +20,69 @@ function ACA(; return ACA(rowpivoting, columnpivoting, convergence) end +function (aca::ACA)(K::AbstractMatrix, ivec::Vector{Int}, jvec::Vector{Int}) + return ACA(aca.rowpivoting(ivec), aca.columnpivoting(jvec), aca.convergence()) +end + +function (aca::ACA{RP,CP,CC})( + K::AbstractMatrix, ivec::Vector{Int}, jvec::Vector{Int} +) where {RP<:PivStrat,CP<:PivStrat,CC<:RandomSampling} + return ACA(aca.rowpivoting(ivec), aca.columnpivoting(jvec), aca.convergence(ivec, jvec)) +end + +function (aca::ACA{RP,CP,CC})( + K::AbstractMatrix, ivec::Vector{Int}, jvec::Vector{Int} +) where {RP<:RandomSamplingPivoting,CP<:PivStrat,CC<:RandomSampling} + convergence = aca.convergence(K, ivec, jvec) + return ACA(aca.rowpivoting(convergence), aca.columnpivoting(jvec), convergence) +end + +function (aca::ACA{RP,CP,CC})( + K::AbstractMatrix, ivec::Vector{Int}, jvec::Vector{Int} +) where {RP<:PivStrat,CP<:RandomSamplingPivoting,CC<:RandomSampling} + convergence = aca.convergence(K, ivec, jvec) + return ACA(aca.rowpivoting(ivec), aca.columnpivoting(convergence), convergence) +end + +function (aca::ACA{RP,CP,CC})( + K::AbstractMatrix, ivec::Vector{Int}, jvec::Vector{Int} +) where {RP<:RandomSamplingPivoting,CP<:RandomSamplingPivoting,CC<:RandomSampling} + convergence = aca.convergence(K, ivec, jvec) + return ACA(aca.rowpivoting(convergence), aca.columnpivoting(convergence), convergence) +end + +function (aca::ACA{RP,CP,CC})( + K::AbstractMatrix, ivec::Vector{Int}, jvec::Vector{Int} +) where {RP<:PivStrat,CP<:CombinedPivStrat,CC<:CombinedConvCrit} + convergence = aca.convergence(K, ivec, jvec) + return ACA(aca.rowpivoting(ivec), aca.columnpivoting(convergence, jvec), convergence) +end + +function (aca::ACA{RP,CP,CC})( + K::AbstractMatrix, ivec::Vector{Int}, jvec::Vector{Int} +) where {RP<:CombinedPivStrat,CP,CC<:CombinedConvCrit} + convergence = aca.convergence(K, ivec, jvec) + return ACA(aca.rowpivoting(convergence, ivec), aca.columnpivoting(jvec), convergence) +end + +function (aca::ACA{RP,CP,CC})( + K::AbstractMatrix, ivec::Vector{Int}, jvec::Vector{Int} +) where {RP<:CombinedPivStrat,CP<:CombinedPivStrat,CC<:CombinedConvCrit} + convergence = aca.convergence(K, ivec, jvec) + return ACA( + aca.rowpivoting(convergence, ivec), + aca.columnpivoting(convergence, jvec), + convergence, + ) +end + function (aca::ACA)( A, rowbuffer::AbstractMatrix{K}, colbuffer::AbstractMatrix{K}, maxrank::Int; - rowidcs=Vector(1:size(colbuffer, 1)), - colidcs=Vector(1:size(rowbuffer, 2)), + rowidcs = Vector(1:size(colbuffer, 1)), + colidcs = Vector(1:size(rowbuffer, 2)), ) where {K} rows = Int[1] cols = Int[] diff --git a/src/convergence/combinedconvcrit.jl b/src/convergence/combinedconvcrit.jl index 0d1f4a1..f1eccac 100644 --- a/src/convergence/combinedconvcrit.jl +++ b/src/convergence/combinedconvcrit.jl @@ -19,3 +19,17 @@ function (convcrit::CombinedConvCrit)( end return npivot, any(convcrit.isconverged) end + +function (convcrit::CombinedConvCrit)( + K::AbstractMatrix, ivec::Vector{Int}, jvec::Vector{Int} +) + curr_crits = Vector{ConvCrit}(undef, length(convcrit.crits)) + for (i, crit) in enumerate(convcrit.crits) + if isa(crit, RandomSampling) + curr_crits[i] = crit(K, ivec, jvec) + else + curr_crits[i] = crit() + end + end + return CombinedConvCrit(curr_crits, zeros(Bool, length(curr_crits))) +end diff --git a/src/convergence/estimation.jl b/src/convergence/estimation.jl index 7285807..b56495e 100644 --- a/src/convergence/estimation.jl +++ b/src/convergence/estimation.jl @@ -7,7 +7,10 @@ function FNormEstimator(tolerance::F) where {F} return FNormEstimator(F(0.0), tolerance) end -(::FNormEstimator{F})() where {F} = FNormEstimator(0.0, 1e-4) +function (cc::FNormEstimator{F})() where {F} + return FNormEstimator(0.0, cc.tol) +end + tolerance(cc::FNormEstimator) = cc.tol function (convcrit::FNormEstimator{F})( diff --git a/src/convergence/extrapolation.jl b/src/convergence/extrapolation.jl index 5721560..fd93f27 100644 --- a/src/convergence/extrapolation.jl +++ b/src/convergence/extrapolation.jl @@ -7,7 +7,9 @@ function FNormExtrapolator(tol::F) where {F} return FNormExtrapolator(F[], FNormEstimator(F(0.0), tol)) end -(::FNormExtrapolator{F})() where {F} = FNormExtrapolator(F[], FNormEstimator(0.0, 1e-4)) +function (cc::FNormExtrapolator{F})() where {F} + return FNormExtrapolator(F[], FNormEstimator(0.0, cc.estimator.tol)) +end tolerance(cc::FNormExtrapolator) = cc.estimator.tol # ACA diff --git a/src/convergence/randomsampling.jl b/src/convergence/randomsampling.jl index f7830a2..856a806 100644 --- a/src/convergence/randomsampling.jl +++ b/src/convergence/randomsampling.jl @@ -14,6 +14,12 @@ function RandomSampling( ) where {K,F<:Real} return RandomSampling{F,K}(F(0.0), nsamples, factor, zeros(Int, 0, 0), zeros(K, 0), tol) end +function (cc::RandomSampling)(K::AbstractMatrix, ivec::Vector{Int}, jvec::Vector{Int}) + nsamples = cc.nsamples == 0 ? cc.factor * (length(ivec) + length(jvec)) : cc.nsamples + indices = hcat(rand(1:length(ivec), nsamples), rand(1:length(jvec), nsamples)) + rest = [K[rc[1], rc[2]][1] for rc in eachrow(indices)] + return RandomSampling(0.0, nsamples, cc.factor, indices, rest, cc.tol) +end tolerance(cc::RandomSampling) = cc.tol diff --git a/src/pivoting/combinedpivstrat.jl b/src/pivoting/combinedpivstrat.jl index c177f61..d48211f 100644 --- a/src/pivoting/combinedpivstrat.jl +++ b/src/pivoting/combinedpivstrat.jl @@ -9,7 +9,7 @@ function (pivstrat::CombinedPivStrat)() end function (pivstrat::CombinedPivStrat)(rc::AbstractArray) - #öprintln(pivstrat.convcrit.isconverged) + #println(pivstrat.convcrit.isconverged) length(pivstrat.strats) > length(pivstrat.convcrit.isconverged) && push!(pivstrat.convcrit.isconverged, false) for (i, conv) in enumerate(pivstrat.convcrit.isconverged) @@ -18,3 +18,15 @@ function (pivstrat::CombinedPivStrat)(rc::AbstractArray) return pivstrat.strats[i](rc) end end + +function (pivstrat::CombinedPivStrat)(convergence::CombinedConvCrit, idcs::Vector{Int}) + curr_strats = Vector{PivStrat}(undef, length(pivstrat.strats)) + for (i, strat) in enumerate(pivstrat.strats) + if isa(strat, RandomSamplingPivoting) + curr_strats[i] = strat(convergence.crits[i]) + else + curr_strats[i] = strat(idcs) + end + end + return CombinedPivStrat(convergence, curr_strats) +end diff --git a/src/pivoting/maxvalue.jl b/src/pivoting/maxvalue.jl index a6ce68c..d477b6a 100644 --- a/src/pivoting/maxvalue.jl +++ b/src/pivoting/maxvalue.jl @@ -6,6 +6,7 @@ end MaximumValue() = MaximumValue(Bool[]) (::MaximumValue)(len::Int) = MaximumValue(zeros(Bool, len)) +(::MaximumValue)(ivec::Vector{Int}) = MaximumValue(zeros(Bool, length(ivec))) function (pivstrat::MaximumValue)() pivstrat.usedidcs[1] = true diff --git a/src/pivoting/randomsampling.jl b/src/pivoting/randomsampling.jl index 8874e43..3228fcd 100644 --- a/src/pivoting/randomsampling.jl +++ b/src/pivoting/randomsampling.jl @@ -6,3 +6,7 @@ end function (pivstrat::RandomSamplingPivoting{F,K})(::AbstractArray) where {F<:Real,K} return pivstrat.convcrit.indices[argmax(pivstrat.convcrit.rest), pivstrat.rc] end + +function (pivstrat::RandomSamplingPivoting{F,K})(convergence::ConvCrit) where {F<:Real,K} + return RandomSamplingPivoting(convergence, pivstrat.rc) +end diff --git a/test/runtests.jl b/test/runtests.jl index f514e67..d6088f5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,6 +5,7 @@ using AdaptiveCrossApproximation include("test_aca.jl") include("test_convergence.jl") include("test_pivoting.jl") + include("test_ACAHMatrices.jl") end @testitem "Code quality (Aqua.jl)" begin @@ -19,7 +20,7 @@ end @testitem "Code formatting (JuliaFormatter.jl)" begin using JuliaFormatter pkgpath = pkgdir(AdaptiveCrossApproximation) - @test JuliaFormatter.format(pkgpath, overwrite=false) + @test JuliaFormatter.format(pkgpath, overwrite=true) end @run_package_tests verbose = true diff --git a/test/test_ACAHMatrices.jl b/test/test_ACAHMatrices.jl index c01bb8c..16ae10f 100644 --- a/test/test_ACAHMatrices.jl +++ b/test/test_ACAHMatrices.jl @@ -1,80 +1,87 @@ -using BEAST using HMatrices -using CompScienceMeshes +using StaticArrays using AdaptiveCrossApproximation using Random using Test using LinearAlgebra -Γ = meshsphere(1.0, 0.2) -op = Helmholtz3D.singlelayer() -spaceX = lagrangecxd0(Γ) -dim = length(spaceX.pos) -xclt = ClusterTree(spaceX.pos) +const Point3D = SVector{3, Float64} +X = rand(Point3D, 1000) +Y = rand(Point3D, 1000) +const k = 2π +function G(buf, irange, jrange) + for (a, i) in enumerate(irange), (b, j) in enumerate(jrange) + buf[a, b] = exp(-k * norm(X[i] - Y[j])) + end + return buf +end + +function G(x, y) + return exp(-k * norm(x - y)) +end -K = HMatrices.KernelMatrix(op, spaceX, spaceX); +K = KernelMatrix(G, X, Y) +dim = size(K, 2) rtols = [1e-2, 1e-4, 1e-6, 1e-8, 1e-10, 1e-12, 1e-14] tst_vec = rand(dim) -fullmat = BEAST.assemble(op, spaceX, spaceX) -trueResult = fullmat * tst_vec; +trueResult = K * tst_vec; # PivStrat :: MaximumValue for rtol in rtols - aca = ACA(; convergence=FNormEstimator(0.0, rtol)) - hmat = HMatrices.assemble_hmatrix(K; comp=aca) - @test size(hmat, 1) == size(fullmat, 1) - @test size(hmat, 2) == size(fullmat, 2) - @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol + aca = ACA(; convergence = FNormEstimator(0.0, rtol)) + hmat = HMatrices.assemble_hmatrix(K; comp = aca) + @test size(hmat, 1) == size(K, 1) + @test size(hmat, 2) == size(K, 2) + @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol end # PivStrat :: leja2 for rtol in rtols - aca = ACA(; rowpivoting=Leja2(spaceX.pos), convergence=FNormEstimator(0.0, rtol)) - hmat = HMatrices.assemble_hmatrix(K; comp=aca) - @test size(hmat, 1) == size(fullmat, 1) - @test size(hmat, 2) == size(fullmat, 2) - @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol + aca = ACA(; rowpivoting = Leja2(X), convergence = FNormEstimator(0.0, rtol)) + hmat = HMatrices.assemble_hmatrix(K; comp = aca) + @test size(hmat, 1) == size(K, 1) + @test size(hmat, 2) == size(K, 2) + @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol - aca = ACA(; columnpivoting=Leja2(spaceX.pos), convergence=FNormEstimator(0.0, rtol)) - hmat = HMatrices.assemble_hmatrix(K; comp=aca) - @test size(hmat, 1) == size(fullmat, 1) - @test size(hmat, 2) == size(fullmat, 2) - @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol + aca = ACA(; columnpivoting = Leja2(Y), convergence = FNormEstimator(0.0, rtol)) + hmat = HMatrices.assemble_hmatrix(K; comp = aca) + @test size(hmat, 1) == size(K, 1) + @test size(hmat, 2) == size(K, 2) + @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol end # PivStrat :: FillDistance for rtol in rtols - aca = ACA(; rowpivoting=FillDistance(spaceX.pos), convergence=FNormEstimator(0.0, rtol)) - hmat = HMatrices.assemble_hmatrix(K; comp=aca) - @test size(hmat, 1) == size(fullmat, 1) - @test size(hmat, 2) == size(fullmat, 2) - @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol + aca = ACA(; rowpivoting = FillDistance(X), convergence = FNormEstimator(0.0, rtol)) + hmat = HMatrices.assemble_hmatrix(K; comp = aca) + @test size(hmat, 1) == size(K, 1) + @test size(hmat, 2) == size(K, 2) + @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol - aca = ACA(; - columnpivoting=FillDistance(spaceX.pos), convergence=FNormEstimator(0.0, rtol) - ) - hmat = HMatrices.assemble_hmatrix(K; comp=aca) - @test size(hmat, 1) == size(fullmat, 1) - @test size(hmat, 2) == size(fullmat, 2) - @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol + aca = ACA(; columnpivoting = FillDistance(Y), convergence = FNormEstimator(0.0, rtol)) + hmat = HMatrices.assemble_hmatrix(K; comp = aca) + @test size(hmat, 1) == size(K, 1) + @test size(hmat, 2) == size(K, 2) + @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol end # PivStrat :: CombinedPivStrat for rtol in [1e-2, 1e-4, 1e-6, 1e-8, 1e-10]#, 1e-12, 1e-14] - Random.seed!(1) - cc1 = AdaptiveCrossApproximation.FNormEstimator(rtol) - cc2 = AdaptiveCrossApproximation.RandomSampling( - Float64; nsamples=100, factor=1.0, tol=rtol - ) - convergence = AdaptiveCrossApproximation.CombinedConvCrit([cc1, cc2], zeros(Bool, 2)) + Random.seed!(1) + cc1 = AdaptiveCrossApproximation.FNormEstimator(rtol) + cc2 = AdaptiveCrossApproximation.RandomSampling( + Float64; nsamples = 100, factor = 1.0, tol = rtol, + ) + convergence = AdaptiveCrossApproximation.CombinedConvCrit([cc1, cc2], zeros(Bool, 2)) - ps1 = MaximumValue() - ps2 = AdaptiveCrossApproximation.RandomSamplingPivoting(cc2, 1) - rp = AdaptiveCrossApproximation.CombinedPivStrat(convergence, [ps1, ps2]) - aca = ACA(; rowpivoting=rp, convergence=convergence) - hmat = HMatrices.assemble_hmatrix(K; comp=aca) - @test size(hmat, 1) == size(fullmat, 1) - @test size(hmat, 2) == size(fullmat, 2) - @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol + ps1 = MaximumValue() + ps2 = AdaptiveCrossApproximation.RandomSamplingPivoting(cc2, 1) + rp = AdaptiveCrossApproximation.CombinedPivStrat(convergence, [ps1, ps2]) + aca = ACA(; rowpivoting = rp, convergence = convergence) + hmat = HMatrices.assemble_hmatrix(K; comp = aca) + @test size(hmat, 1) == size(K, 1) + @test size(hmat, 2) == size(K, 2) + @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol end + From aec185dcb124a7780b2ade817b2788d8de61e7c4 Mon Sep 17 00:00:00 2001 From: aliat Date: Tue, 19 Aug 2025 11:47:46 +0200 Subject: [PATCH 3/4] Pivoting passed tests --- Project.toml | 6 ++- src/aca.jl | 66 +++++++++-------------- src/pivoting/combinedpivstrat.jl | 1 - test/runtests.jl | 2 +- test/test_ACAHMatrices.jl | 90 ++++++++++++++++---------------- 5 files changed, 75 insertions(+), 90 deletions(-) diff --git a/Project.toml b/Project.toml index f8e3b78..6c65bd3 100644 --- a/Project.toml +++ b/Project.toml @@ -4,11 +4,13 @@ authors = ["JoshuaTetzner "] version = "0.1.0" [deps] -HMatrices = "8646bddf-ab1c-4fa7-9c51-ba187d647618" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +[weakdeps] +HMatrices = "8646bddf-ab1c-4fa7-9c51-ba187d647618" + [extensions] ACAHMatrices = ["HMatrices"] @@ -29,4 +31,4 @@ TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a" TestItems = "1c621080-faea-4a02-84b6-bbd5e436b8fe" [targets] -test = ["Aqua", "JET", "JuliaFormatter", "Test", "TestItemRunner", "TestItems", "Random"] +test = ["Aqua", "JET", "JuliaFormatter", "Test", "TestItemRunner", "TestItems", "Random", "HMatrices"] diff --git a/src/aca.jl b/src/aca.jl index 93f810f..8461d28 100644 --- a/src/aca.jl +++ b/src/aca.jl @@ -24,56 +24,40 @@ function (aca::ACA)(K::AbstractMatrix, ivec::Vector{Int}, jvec::Vector{Int}) return ACA(aca.rowpivoting(ivec), aca.columnpivoting(jvec), aca.convergence()) end -function (aca::ACA{RP,CP,CC})( - K::AbstractMatrix, ivec::Vector{Int}, jvec::Vector{Int} -) where {RP<:PivStrat,CP<:PivStrat,CC<:RandomSampling} - return ACA(aca.rowpivoting(ivec), aca.columnpivoting(jvec), aca.convergence(ivec, jvec)) -end - -function (aca::ACA{RP,CP,CC})( +function (aca::ACA{PS,PS,CC})( K::AbstractMatrix, ivec::Vector{Int}, jvec::Vector{Int} -) where {RP<:RandomSamplingPivoting,CP<:PivStrat,CC<:RandomSampling} +) where {PS<:PivStrat,CC<:RandomSampling} convergence = aca.convergence(K, ivec, jvec) - return ACA(aca.rowpivoting(convergence), aca.columnpivoting(jvec), convergence) -end - -function (aca::ACA{RP,CP,CC})( - K::AbstractMatrix, ivec::Vector{Int}, jvec::Vector{Int} -) where {RP<:PivStrat,CP<:RandomSamplingPivoting,CC<:RandomSampling} - convergence = aca.convergence(K, ivec, jvec) - return ACA(aca.rowpivoting(ivec), aca.columnpivoting(convergence), convergence) -end - -function (aca::ACA{RP,CP,CC})( - K::AbstractMatrix, ivec::Vector{Int}, jvec::Vector{Int} -) where {RP<:RandomSamplingPivoting,CP<:RandomSamplingPivoting,CC<:RandomSampling} - convergence = aca.convergence(K, ivec, jvec) - return ACA(aca.rowpivoting(convergence), aca.columnpivoting(convergence), convergence) -end + rowpivoting = if aca.rowpivoting isa RandomSamplingPivoting + aca.rowpivoting(convergence) + else + aca.rowpivoting(ivec) + end + columnpivoting = if aca.columnpivoting isa RandomSamplingPivoting + aca.columnpivoting(convergence) + else + aca.columnpivoting(jvec) + end -function (aca::ACA{RP,CP,CC})( - K::AbstractMatrix, ivec::Vector{Int}, jvec::Vector{Int} -) where {RP<:PivStrat,CP<:CombinedPivStrat,CC<:CombinedConvCrit} - convergence = aca.convergence(K, ivec, jvec) - return ACA(aca.rowpivoting(ivec), aca.columnpivoting(convergence, jvec), convergence) + return ACA(rowpivoting, columnpivoting, convergence) end function (aca::ACA{RP,CP,CC})( K::AbstractMatrix, ivec::Vector{Int}, jvec::Vector{Int} -) where {RP<:CombinedPivStrat,CP,CC<:CombinedConvCrit} +) where {RP,CP,CC<:CombinedConvCrit} convergence = aca.convergence(K, ivec, jvec) - return ACA(aca.rowpivoting(convergence, ivec), aca.columnpivoting(jvec), convergence) -end + rowpivoting = if aca.rowpivoting isa CombinedPivStrat + aca.rowpivoting(convergence, ivec) + else + aca.rowpivoting(ivec) + end + columnpivoting = if aca.columnpivoting isa CombinedPivStrat + aca.columnpivoting(convergence, jvec) + else + aca.columnpivoting(jvec) + end -function (aca::ACA{RP,CP,CC})( - K::AbstractMatrix, ivec::Vector{Int}, jvec::Vector{Int} -) where {RP<:CombinedPivStrat,CP<:CombinedPivStrat,CC<:CombinedConvCrit} - convergence = aca.convergence(K, ivec, jvec) - return ACA( - aca.rowpivoting(convergence, ivec), - aca.columnpivoting(convergence, jvec), - convergence, - ) + return ACA(rowpivoting, columnpivoting, convergence) end function (aca::ACA)( diff --git a/src/pivoting/combinedpivstrat.jl b/src/pivoting/combinedpivstrat.jl index d48211f..95fa6df 100644 --- a/src/pivoting/combinedpivstrat.jl +++ b/src/pivoting/combinedpivstrat.jl @@ -9,7 +9,6 @@ function (pivstrat::CombinedPivStrat)() end function (pivstrat::CombinedPivStrat)(rc::AbstractArray) - #println(pivstrat.convcrit.isconverged) length(pivstrat.strats) > length(pivstrat.convcrit.isconverged) && push!(pivstrat.convcrit.isconverged, false) for (i, conv) in enumerate(pivstrat.convcrit.isconverged) diff --git a/test/runtests.jl b/test/runtests.jl index d6088f5..205411f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -20,7 +20,7 @@ end @testitem "Code formatting (JuliaFormatter.jl)" begin using JuliaFormatter pkgpath = pkgdir(AdaptiveCrossApproximation) - @test JuliaFormatter.format(pkgpath, overwrite=true) + @test JuliaFormatter.format(pkgpath, overwrite=false) end @run_package_tests verbose = true diff --git a/test/test_ACAHMatrices.jl b/test/test_ACAHMatrices.jl index 16ae10f..83cc651 100644 --- a/test/test_ACAHMatrices.jl +++ b/test/test_ACAHMatrices.jl @@ -5,19 +5,19 @@ using Random using Test using LinearAlgebra -const Point3D = SVector{3, Float64} +const Point3D = SVector{3,Float64} X = rand(Point3D, 1000) Y = rand(Point3D, 1000) const k = 2π function G(buf, irange, jrange) - for (a, i) in enumerate(irange), (b, j) in enumerate(jrange) - buf[a, b] = exp(-k * norm(X[i] - Y[j])) - end - return buf + for (a, i) in enumerate(irange), (b, j) in enumerate(jrange) + buf[a, b] = exp(-k * norm(X[i] - Y[j])) + end + return buf end function G(x, y) - return exp(-k * norm(x - y)) + return exp(-k * norm(x - y)) end K = KernelMatrix(G, X, Y) @@ -30,58 +30,58 @@ trueResult = K * tst_vec; # PivStrat :: MaximumValue for rtol in rtols - aca = ACA(; convergence = FNormEstimator(0.0, rtol)) - hmat = HMatrices.assemble_hmatrix(K; comp = aca) - @test size(hmat, 1) == size(K, 1) - @test size(hmat, 2) == size(K, 2) - @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol + aca = ACA(; convergence=FNormEstimator(0.0, rtol)) + hmat = HMatrices.assemble_hmatrix(K; comp=aca) + @test size(hmat, 1) == size(K, 1) + @test size(hmat, 2) == size(K, 2) + @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol end # PivStrat :: leja2 for rtol in rtols - aca = ACA(; rowpivoting = Leja2(X), convergence = FNormEstimator(0.0, rtol)) - hmat = HMatrices.assemble_hmatrix(K; comp = aca) - @test size(hmat, 1) == size(K, 1) - @test size(hmat, 2) == size(K, 2) - @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol + aca = ACA(; rowpivoting=Leja2(X), convergence=FNormEstimator(0.0, rtol)) + hmat = HMatrices.assemble_hmatrix(K; comp=aca) + @test size(hmat, 1) == size(K, 1) + @test size(hmat, 2) == size(K, 2) + @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol - aca = ACA(; columnpivoting = Leja2(Y), convergence = FNormEstimator(0.0, rtol)) - hmat = HMatrices.assemble_hmatrix(K; comp = aca) - @test size(hmat, 1) == size(K, 1) - @test size(hmat, 2) == size(K, 2) - @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol + aca = ACA(; columnpivoting=Leja2(Y), convergence=FNormEstimator(0.0, rtol)) + hmat = HMatrices.assemble_hmatrix(K; comp=aca) + @test size(hmat, 1) == size(K, 1) + @test size(hmat, 2) == size(K, 2) + @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol end # PivStrat :: FillDistance for rtol in rtols - aca = ACA(; rowpivoting = FillDistance(X), convergence = FNormEstimator(0.0, rtol)) - hmat = HMatrices.assemble_hmatrix(K; comp = aca) - @test size(hmat, 1) == size(K, 1) - @test size(hmat, 2) == size(K, 2) - @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol + aca = ACA(; rowpivoting=FillDistance(X), convergence=FNormEstimator(0.0, rtol)) + hmat = HMatrices.assemble_hmatrix(K; comp=aca) + @test size(hmat, 1) == size(K, 1) + @test size(hmat, 2) == size(K, 2) + @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol - aca = ACA(; columnpivoting = FillDistance(Y), convergence = FNormEstimator(0.0, rtol)) - hmat = HMatrices.assemble_hmatrix(K; comp = aca) - @test size(hmat, 1) == size(K, 1) - @test size(hmat, 2) == size(K, 2) - @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol + aca = ACA(; columnpivoting=FillDistance(Y), convergence=FNormEstimator(0.0, rtol)) + hmat = HMatrices.assemble_hmatrix(K; comp=aca) + @test size(hmat, 1) == size(K, 1) + @test size(hmat, 2) == size(K, 2) + @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol end # PivStrat :: CombinedPivStrat for rtol in [1e-2, 1e-4, 1e-6, 1e-8, 1e-10]#, 1e-12, 1e-14] - Random.seed!(1) - cc1 = AdaptiveCrossApproximation.FNormEstimator(rtol) - cc2 = AdaptiveCrossApproximation.RandomSampling( - Float64; nsamples = 100, factor = 1.0, tol = rtol, - ) - convergence = AdaptiveCrossApproximation.CombinedConvCrit([cc1, cc2], zeros(Bool, 2)) + Random.seed!(1) + cc1 = AdaptiveCrossApproximation.FNormEstimator(rtol) + cc2 = AdaptiveCrossApproximation.RandomSampling( + Float64; nsamples=100, factor=1.0, tol=rtol + ) + convergence = AdaptiveCrossApproximation.CombinedConvCrit([cc1, cc2], zeros(Bool, 2)) - ps1 = MaximumValue() - ps2 = AdaptiveCrossApproximation.RandomSamplingPivoting(cc2, 1) - rp = AdaptiveCrossApproximation.CombinedPivStrat(convergence, [ps1, ps2]) - aca = ACA(; rowpivoting = rp, convergence = convergence) - hmat = HMatrices.assemble_hmatrix(K; comp = aca) - @test size(hmat, 1) == size(K, 1) - @test size(hmat, 2) == size(K, 2) - @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol + ps1 = MaximumValue() + ps2 = AdaptiveCrossApproximation.RandomSamplingPivoting(cc2, 1) + rp = AdaptiveCrossApproximation.CombinedPivStrat(convergence, [ps1, ps2]) + aca = ACA(; rowpivoting=rp, convergence=convergence) + hmat = HMatrices.assemble_hmatrix(K; comp=aca) + @test size(hmat, 1) == size(K, 1) + @test size(hmat, 2) == size(K, 2) + @test norm(hmat * tst_vec - trueResult) / norm(trueResult) ≈ 0 atol = rtol end From b12634f7354a8a32689f66975f584cf057623142 Mon Sep 17 00:00:00 2001 From: aliat Date: Tue, 19 Aug 2025 21:42:13 +0200 Subject: [PATCH 4/4] using Abstractarray(Int) --- ext/ACAHMatrices/hmatrices_interface.jl | 2 +- src/aca.jl | 26 +++++++++++++------------ src/convergence/combinedconvcrit.jl | 4 ++-- src/convergence/randomsampling.jl | 10 +++++++--- src/pivoting/combinedpivstrat.jl | 4 +++- src/pivoting/filldistance.jl | 2 +- src/pivoting/lejapoints.jl | 2 +- src/pivoting/maxvalue.jl | 2 +- 8 files changed, 30 insertions(+), 22 deletions(-) diff --git a/ext/ACAHMatrices/hmatrices_interface.jl b/ext/ACAHMatrices/hmatrices_interface.jl index f10e659..c08a9ac 100644 --- a/ext/ACAHMatrices/hmatrices_interface.jl +++ b/ext/ACAHMatrices/hmatrices_interface.jl @@ -19,7 +19,7 @@ function (aca::ACA)( bufs.B.data[1:(length(irange) * maxrank)], (length(irange), maxrank) ) - aca = aca(K, Vector(irange), Vector(jrange)) + aca = aca(K, irange, jrange) npivots, U, V = aca(K, rowBuffer, colBuffer, maxrank; rowidcs=irange, colidcs=jrange) return HMatrices.RkMatrix( colBuffer[:, 1:npivots], Matrix(transpose(rowBuffer[1:npivots, :])) diff --git a/src/aca.jl b/src/aca.jl index 8461d28..1a90c30 100644 --- a/src/aca.jl +++ b/src/aca.jl @@ -20,41 +20,43 @@ function ACA(; return ACA(rowpivoting, columnpivoting, convergence) end -function (aca::ACA)(K::AbstractMatrix, ivec::Vector{Int}, jvec::Vector{Int}) - return ACA(aca.rowpivoting(ivec), aca.columnpivoting(jvec), aca.convergence()) +function (aca::ACA)( + K::AbstractMatrix, rowidcs::AbstractArray{Int}, colidcs::AbstractArray{Int} +) + return ACA(aca.rowpivoting(rowidcs), aca.columnpivoting(colidcs), aca.convergence()) end function (aca::ACA{PS,PS,CC})( - K::AbstractMatrix, ivec::Vector{Int}, jvec::Vector{Int} + K::AbstractMatrix, rowidcs::AbstractArray{Int}, colidcs::AbstractArray{Int} ) where {PS<:PivStrat,CC<:RandomSampling} - convergence = aca.convergence(K, ivec, jvec) + convergence = aca.convergence(K, rowidcs, colidcs) rowpivoting = if aca.rowpivoting isa RandomSamplingPivoting aca.rowpivoting(convergence) else - aca.rowpivoting(ivec) + aca.rowpivoting(rowidcs) end columnpivoting = if aca.columnpivoting isa RandomSamplingPivoting aca.columnpivoting(convergence) else - aca.columnpivoting(jvec) + aca.columnpivoting(colidcs) end return ACA(rowpivoting, columnpivoting, convergence) end function (aca::ACA{RP,CP,CC})( - K::AbstractMatrix, ivec::Vector{Int}, jvec::Vector{Int} + K::AbstractMatrix, rowidcs::AbstractArray{Int}, colidcs::AbstractArray{Int} ) where {RP,CP,CC<:CombinedConvCrit} - convergence = aca.convergence(K, ivec, jvec) + convergence = aca.convergence(K, rowidcs, colidcs) rowpivoting = if aca.rowpivoting isa CombinedPivStrat - aca.rowpivoting(convergence, ivec) + aca.rowpivoting(convergence, rowidcs) else - aca.rowpivoting(ivec) + aca.rowpivoting(rowidcs) end columnpivoting = if aca.columnpivoting isa CombinedPivStrat - aca.columnpivoting(convergence, jvec) + aca.columnpivoting(convergence, colidcs) else - aca.columnpivoting(jvec) + aca.columnpivoting(colidcs) end return ACA(rowpivoting, columnpivoting, convergence) diff --git a/src/convergence/combinedconvcrit.jl b/src/convergence/combinedconvcrit.jl index f1eccac..095395b 100644 --- a/src/convergence/combinedconvcrit.jl +++ b/src/convergence/combinedconvcrit.jl @@ -21,12 +21,12 @@ function (convcrit::CombinedConvCrit)( end function (convcrit::CombinedConvCrit)( - K::AbstractMatrix, ivec::Vector{Int}, jvec::Vector{Int} + K::AbstractMatrix, rowidcs::AbstractArray{Int}, colidcs::AbstractArray{Int} ) curr_crits = Vector{ConvCrit}(undef, length(convcrit.crits)) for (i, crit) in enumerate(convcrit.crits) if isa(crit, RandomSampling) - curr_crits[i] = crit(K, ivec, jvec) + curr_crits[i] = crit(K, rowidcs, colidcs) else curr_crits[i] = crit() end diff --git a/src/convergence/randomsampling.jl b/src/convergence/randomsampling.jl index 856a806..6b4cc8d 100644 --- a/src/convergence/randomsampling.jl +++ b/src/convergence/randomsampling.jl @@ -14,9 +14,13 @@ function RandomSampling( ) where {K,F<:Real} return RandomSampling{F,K}(F(0.0), nsamples, factor, zeros(Int, 0, 0), zeros(K, 0), tol) end -function (cc::RandomSampling)(K::AbstractMatrix, ivec::Vector{Int}, jvec::Vector{Int}) - nsamples = cc.nsamples == 0 ? cc.factor * (length(ivec) + length(jvec)) : cc.nsamples - indices = hcat(rand(1:length(ivec), nsamples), rand(1:length(jvec), nsamples)) +function (cc::RandomSampling)( + K::AbstractMatrix, rowidcs::AbstractArray{Int}, colidcs::AbstractArray{Int} +) + rowlen = length(rowidcs) + collen = length(colidcs) + nsamples = cc.nsamples == 0 ? cc.factor * (rowlen + collen) : cc.nsamples + indices = hcat(rand(1:rowlen, nsamples), rand(1:collen, nsamples)) rest = [K[rc[1], rc[2]][1] for rc in eachrow(indices)] return RandomSampling(0.0, nsamples, cc.factor, indices, rest, cc.tol) end diff --git a/src/pivoting/combinedpivstrat.jl b/src/pivoting/combinedpivstrat.jl index 95fa6df..7ecc1a3 100644 --- a/src/pivoting/combinedpivstrat.jl +++ b/src/pivoting/combinedpivstrat.jl @@ -18,7 +18,9 @@ function (pivstrat::CombinedPivStrat)(rc::AbstractArray) end end -function (pivstrat::CombinedPivStrat)(convergence::CombinedConvCrit, idcs::Vector{Int}) +function (pivstrat::CombinedPivStrat)( + convergence::CombinedConvCrit, idcs::AbstractArray{Int} +) curr_strats = Vector{PivStrat}(undef, length(pivstrat.strats)) for (i, strat) in enumerate(pivstrat.strats) if isa(strat, RandomSamplingPivoting) diff --git a/src/pivoting/filldistance.jl b/src/pivoting/filldistance.jl index 25165a1..58f3bc1 100644 --- a/src/pivoting/filldistance.jl +++ b/src/pivoting/filldistance.jl @@ -8,7 +8,7 @@ function FillDistance(pos::Vector{SVector{D,F}}) where {D,F<:Real} return FillDistance{D,F}(F[], pos) end -function (pivstrat::FillDistance{D,F})(idcs::Vector{Int}) where {D,F} +function (pivstrat::FillDistance{D,F})(idcs::AbstractArray{Int}) where {D,F} return FillDistance{D,F}(zeros(F, length(idcs)), pivstrat.pos[idcs]) end diff --git a/src/pivoting/lejapoints.jl b/src/pivoting/lejapoints.jl index 697f92e..dcdd885 100644 --- a/src/pivoting/lejapoints.jl +++ b/src/pivoting/lejapoints.jl @@ -16,7 +16,7 @@ function Leja2(pos::Vector{SVector{D,F}}) where {D,F} return Leja2{D,F}(F[], pos) end -function (pivstrat::Leja2{D,F})(idcs::Vector{Int}) where {D,F} +function (pivstrat::Leja2{D,F})(idcs::AbstractArray{Int}) where {D,F} return Leja2{D,F}(zeros(F, length(idcs)), pivstrat.pos[idcs]) end diff --git a/src/pivoting/maxvalue.jl b/src/pivoting/maxvalue.jl index d477b6a..85205c1 100644 --- a/src/pivoting/maxvalue.jl +++ b/src/pivoting/maxvalue.jl @@ -6,7 +6,7 @@ end MaximumValue() = MaximumValue(Bool[]) (::MaximumValue)(len::Int) = MaximumValue(zeros(Bool, len)) -(::MaximumValue)(ivec::Vector{Int}) = MaximumValue(zeros(Bool, length(ivec))) +(::MaximumValue)(idcs::AbstractArray{Int}) = MaximumValue(zeros(Bool, length(idcs))) function (pivstrat::MaximumValue)() pivstrat.usedidcs[1] = true