fixed a bug in ulvdivide.jl regarding unbalanced cluster trees

Project.toml

@@ -1,7 +1,7 @@
 name = "HssMatrices"
 uuid = "0c1fde87-e0ea-49d6-ac2f-89b694f639e1"
 authors = ["Boris Bonev <[email protected]>"]
-version = "0.1.3"
+version = "0.1.4"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"

src/ulvdivide.jl

@@ -191,8 +191,8 @@ function _unpackadd_rows!(hssA::HssMatrix{T}, hssB::HssMatrix{T}, ktree::BinaryN
 
       hssA.A11 = _unpackadd_rows!(hssA.A11, hssB.A11, ktree.left)
       hssA.A22 = _unpackadd_rows!(hssA.A22, hssB.A22, ktree.right)
-      hssA.sz1 == size(hssA.A11) || error("Didn't expect dimensions to change")
-      hssA.sz2 == size(hssA.A22) || error("Didn't expect dimensions to change")
+      hssA.sz1 == size(hssA.A11) || error("dimensions got screwed up!")
+      hssA.sz2 == size(hssA.A22) || error("dimensions got screwed up!")
     else
       isbranch(ktree) || throw(ArgumentError("didn't expect ktree to be a leaf node"))
       k1 = ktree.left.data; k2 = ktree.right.data
@@ -218,7 +218,17 @@ function _unpackadd_rows!(hssA::HssMatrix{T}, hssB::HssMatrix{T}, ktree::BinaryN
     end
     return hssA
   else
-    error("Expected hssA to be a branch node")
+    if isbranch(hssB)
+      error("Something went wrong. Expected hssB to also be a leaf node")
+    else
+      k, rk = size(hssB.U);
+      k <= size(hssA,1) || throw(DimensionMismatch("expected hssB to have less rows than hssA."))
+      hssA.U = [hssA.U zeros(size(hssA.U,1), rk)]
+      hssA.U[end-k+1:end, end-rk+1:end] .= hssB.U
+      hssA.V = [hssA.V hssB.V]
+      hssA.D[end-k+1:end, :] .= hssA.D[end-k+1:end, :] .+ hssB.D
+      return hssA
+    end
   end
 end
 

test/examples.jl

@@ -8,6 +8,7 @@ using Plots
 K(x,y) = (x-y) > 0 ? 0.001/(x-y) : 2.
 #K(x,y) = (x-y) != 0 ? 1/(x-y) : 1.
 A = [ K(x,y) for x=-1:0.001:1, y=-1:0.001:1];
+C = inv(A);
 
 # test the simple implementation of cluster trees
 m, n = size(A)
@@ -50,16 +51,24 @@ Id(i,j) = Matrix{Float64}(i.*ones(length(j))' .== ones(length(i)).*j')
 IdOp = LinearMap{Float64}(n, n, (y,_,x) -> x, (y,_,x) -> x, (i,j) -> Id(i,j))
 hssI = randcompress(IdOp, ccl, ccl, 0)
 hssC = ldiv!(hssA, hssI)
-println("rel. error in the solution of AX = I: ", norm(full(hssC) - inv(A)) / norm(inv(A)) )
-println("abs. error in the solution of AX = I: ", norm(full(hssC) - inv(A)) )
+println("rel. error in the solution of AX = I: ", norm(full(hssC) - C) / norm(C) )
+println("abs. error in the solution of AX = I: ", norm(full(hssC) - C) )
 
 # test left HSS division
 Id(i,j) = Matrix{Float64}(i.*ones(length(j))' .== ones(length(i)).*j')
 IdOp = LinearMap{Float64}(n, n, (y,_,x) -> x, (y,_,x) -> x, (i,j) -> Id(i,j))
 hssI = randcompress(IdOp, ccl, ccl, 0)
 hssC = rdiv!(hssI, hssA)
-println("rel. error in the solution of XA = I: ", norm(full(hssC) - inv(A)) / norm(inv(A)) )
-println("abs. error in the solution of XA = I: ", norm(full(hssC) - inv(A)) )
+println("rel. error in the solution of XA = I: ", norm(full(hssC) - C) / norm(C) )
+println("abs. error in the solution of XA = I: ", norm(full(hssC) - C) )
+
+# desequilibrate cluster trees and try again
+hssA.A11 = prune_leaves!(hssA.A11)
+hssI = randcompress(IdOp, ccl, ccl, 0)
+hssI.A11 = prune_leaves!(hssI.A11)
+hssC = ldiv!(hssA, hssI)
+println("rel. error in the solution of XA = I: ", norm(full(hssC) - C) / norm(C) )
+println("abs. error in the solution of XA = I: ", norm(full(hssC) - C) )
 
 
 # # test computation of generators
@@ -79,4 +88,4 @@ println("abs. error in the solution of XA = I: ", norm(full(hssC) - inv(A)) )
 # @time x = ulvfactsolve(hssA, b);
 
 # # test plotting
-# plot = plotranks(hssA)
+# plot = plotranks(hssA)

test/runtests.jl

@@ -53,4 +53,8 @@ end;
     Ainv = inv(A)
     @test norm(Ainv - full(hssA\hssI))/norm(Ainv) ≤ c*HssOptions().rtol || norm(Ainv - full(hssA\hssI)) ≤ c*HssOptions().atol
     @test norm(Ainv - full(hssI/hssA))/norm(Ainv) ≤ c*HssOptions().rtol || norm(Ainv - full(hssI/hssA)) ≤ c*HssOptions().atol
+    hssA.A11 = prune_leaves!(hssA.A11)
+    hssI.A11 = prune_leaves!(hssI.A11)
+    @test norm(Ainv - full(hssA\hssI))/norm(Ainv) ≤ c*HssOptions().rtol || norm(Ainv - full(hssA\hssI)) ≤ c*HssOptions().atol
+    @test norm(Ainv - full(hssI/hssA))/norm(Ainv) ≤ c*HssOptions().rtol || norm(Ainv - full(hssI/hssA)) ≤ c*HssOptions().atol
 end