From 06a4f25cc9dbe40d0dbdac0d6e2d3d1f1b50ed2b Mon Sep 17 00:00:00 2001
From: Lyndon White <lyndon.white@invenialabs.co.uk>
Date: Thu, 21 Mar 2019 20:21:28 +0000
Subject: [PATCH 1/2] add large matrix inversion

---
 src/GenericLinearAlgebra.jl |  1 +
 src/inv.jl                  | 68 +++++++++++++++++++++++++++++++++++++
 test/inv.jl                 | 22 ++++++++++++
 test/runtests.jl            |  1 +
 4 files changed, 92 insertions(+)
 create mode 100644 src/inv.jl
 create mode 100644 test/inv.jl

diff --git a/src/GenericLinearAlgebra.jl b/src/GenericLinearAlgebra.jl
index 2ab778b..dea8121 100644
--- a/src/GenericLinearAlgebra.jl
+++ b/src/GenericLinearAlgebra.jl
@@ -12,4 +12,5 @@ include("eigenGeneral.jl")
 include("tridiag.jl")
 include("svd.jl")
 include("rectfullpacked.jl")
+include("inv.jl")
 end
diff --git a/src/inv.jl b/src/inv.jl
new file mode 100644
index 0000000..6b2e495
--- /dev/null
+++ b/src/inv.jl
@@ -0,0 +1,68 @@
+export large_inv
+
+function block_inv(A, B, C, D_inv)
+    size(A, 1) != size(A, 2) && throw(DimensionMismatch("block A is not square."))
+    
+    B_D_inv = B * D_inv
+    pre_inv = inv(A - B_D_inv * C);
+    pre_inv_B_D_inv = pre_inv * B_D_inv
+    D_inv_C = D_inv * C
+    B3 = -D_inv_C * pre_inv
+    B4 = D_inv + D_inv_C * pre_inv_B_D_inv
+        
+    mat_inv = [
+        pre_inv    -pre_inv_B_D_inv
+        B3         B4
+    ]
+    return mat_inv
+end
+
+
+@views function partition_large_mat(mat; max_block_size)
+    size(mat, 1) != size(mat, 2) && throw(DimensionMismatch("Matrix is not square."))
+    Bl_1 = mat[1:max_block_size, 1:max_block_size];
+    Bl_2 = mat[1:max_block_size, max_block_size+1:end];
+    Bl_3 = mat[max_block_size+1:end, 1:max_block_size];
+    Bl_4 = mat[max_block_size+1:end, max_block_size+1:end];
+    return Bl_1, Bl_2, Bl_3, Bl_4
+end
+
+
+function blocks_large_mat(mat::T; max_block) where T<:AbstractMatrix{F} where F
+    if size(mat, 1) <= max_block
+        throw(ArgumentError("The matrix size is smaller than the specified maximum block size."))
+    end
+    # SubMat is the type that `partition_large_mat` returns
+    SubMat = SubArray{F,2,T,Tuple{UnitRange{Int64},UnitRange{Int64}},false}
+    mat_blocks = Tuple{SubMat,SubMat,SubMat,SubMat}[]
+    D = mat
+    while true
+        A, B, C, D = partition_large_mat(D, max_block_size = max_block)
+        push!(mat_blocks, (A, B, C, D))
+        size(D, 1) <= max_block && break
+    end
+    return mat_blocks
+end
+
+
+"""
+    large_inv(mat; max_block_size)
+
+Inverts a large matrix via recursive application of the matrix inversion lemma.
+The matrix is recursively partitioned into blocks until those blocks are below `max_block_size`.
+For large matrices it is significantly faster than `inv`.
+However, it accumulates floating point error faster.
+Particularly when `max_block_size` is small relative to the size of `mat`.
+"""
+function large_inv(mat; max_block_size)
+    blocks = blocks_large_mat(mat, max_block = max_block_size);
+    (A,B,C,D) = pop!(blocks)
+    inverted_mat = block_inv(A, B, C, inv(D));
+
+    @debug "iteration $iter is compeleted."
+    while(!isempty(blocks))
+        (A,B,C,D) = pop!(blocks)
+        inverted_mat = block_inv(A, B, C, inverted_mat);
+    end
+    return inverted_mat
+end
diff --git a/test/inv.jl b/test/inv.jl
new file mode 100644
index 0000000..d3777fa
--- /dev/null
+++ b/test/inv.jl
@@ -0,0 +1,22 @@
+using Test
+using LinearAlgebra
+using GenericLinearAlgebra
+
+@testset "large matrix inv  " begin
+    @testset "Problem dimension ($m,$m) with block size $bz" for
+        (m, bz) in ( # Standard
+                    ( 100,  50), ( 1000,  50), (1000,  500), (10_000, 5_000),
+                     # Nondivisable by block size
+                    ( 100,  64), ( 1000,  151), (1000,  331), (10_000, 6331),
+                   )
+
+        A = rand(m, m)
+        A_inv = large_inv(A; max_block_size=bz)
+        @test A_inv*A - I ≈ zeros(m,m) atol=1e-5
+    end
+
+    @testset "Error paths" begin
+        @test_throws DimensionMismatch large_inv(ones(700,300); max_block_size=100)
+        @test_throws ArgumentError large_inv(ones(100,100); max_block_size=1000)
+    end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 3f10df1..0af0a7a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -9,4 +9,5 @@ using Test
     include("svd.jl")
     include("rectfullpacked.jl")
     include("lapack.jl")
+    include("inv.jl")
 # end

From c5e5a592916453c2b61725a91c967a1f84b31a77 Mon Sep 17 00:00:00 2001
From: Lyndon White <oxinabox@ucc.asn.au>
Date: Mon, 1 Apr 2019 12:17:45 +0100
Subject: [PATCH 2/2] increase tollerance

---
 test/inv.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/inv.jl b/test/inv.jl
index d3777fa..f2a7d6c 100644
--- a/test/inv.jl
+++ b/test/inv.jl
@@ -12,7 +12,7 @@ using GenericLinearAlgebra
 
         A = rand(m, m)
         A_inv = large_inv(A; max_block_size=bz)
-        @test A_inv*A - I ≈ zeros(m,m) atol=1e-5
+        @test A_inv*A - I ≈ zeros(m,m) atol=1e-4
     end
 
     @testset "Error paths" begin