From 787445f37c1606034545b41cd66782188ebee928 Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Fri, 2 Oct 2020 08:25:13 -0400
Subject: [PATCH] mjo/ldlt.py: refactor the one-by-one block-LDLT pivot.

---
 mjo/ldlt.py | 76 +++++++++++++++++++++--------------------------------
 1 file changed, 30 insertions(+), 46 deletions(-)

diff --git a/mjo/ldlt.py b/mjo/ldlt.py
index 5c1e65d..c6c4be3 100644
--- a/mjo/ldlt.py
+++ b/mjo/ldlt.py
@@ -260,8 +260,10 @@ def block_ldlt_naive(A, check_hermitian=False):
         #   [ 1 0 ]
         #   [ 0 B ]
         #
+        # We could still do a pivot_one_by_one() here, but it would
+        # pointlessly subract a bunch of zeros and multiply by one.
         B = A1[1:,1:]
-        P2, L2, D2 = ldlt_naive(B)
+        P2, L2, D2 = block_ldlt_naive(B)
         P1 = block_matrix(2,2, [[ZZ(1), ZZ(0)],
                                 [ZZ(0), P2]])
         L1 = block_matrix(2,2, [[ZZ(1), ZZ(0)],
@@ -270,24 +272,38 @@ def block_ldlt_naive(A, check_hermitian=False):
                                 [ZZ(0), D2]])
         return (P1,L1,D1)
 
-    if A1[0,0].abs() > alpha*omega_1:
-        # Higham step (1)
-        # The top-left entry is our 1x1 pivot.
-        C = A1[1:n,0]
-        B = A1[1:,1:]
+    def pivot_one_by_one(M, c=None):
+        # Perform a one-by-one pivot on "M," swapping row/columns "c".
+        # If "c" is None, no swap is performed.
+        if c is not None:
+            P1 = copy(M.matrix_space().identity_matrix())
+            P1.swap_rows(0,c)
+            M = P1.T * M * P1
 
-        P2, L2, D2 = block_ldlt_naive(B - (C*C.transpose())/A1[0,0])
+        # The top-left entry is now our 1x1 pivot.
+        C = M[1:n,0]
+        B = M[1:,1:]
+
+        P2, L2, D2 = block_ldlt_naive(B - (C*C.transpose())/M[0,0])
+
+        if c is None:
+            P1 = block_matrix(2,2, [[ZZ(1), ZZ(0)],
+                                    [0*C,   P2]])
+        else:
+            P1 = P1*block_matrix(2,2, [[ZZ(1), ZZ(0)],
+                                       [0*C,   P2]])
 
-        P1 = block_matrix(2,2, [[ZZ(1), ZZ(0)],
-                                [0*C,  P2]])
         L1 = block_matrix(2,2, [[ZZ(1),                    ZZ(0)],
-                                [P2.transpose()*C/A1[0,0], L2]])
-        D1 = block_matrix(2,2, [[A1[0,0], ZZ(0)],
-                                [0*C,   D2]])
+                                [P2.transpose()*C/M[0,0], L2]])
+        D1 = block_matrix(2,2, [[M[0,0], ZZ(0)],
+                                [0*C,     D2]])
 
         return (P1,L1,D1)
 
 
+    if A1[0,0].abs() > alpha*omega_1:
+        return pivot_one_by_one(A1)
+
     r = 1 + column_1_subdiag.index(omega_1)
 
     # If the matrix is Hermitian, we need only look at the above-
@@ -295,44 +311,12 @@ def block_ldlt_naive(A, check_hermitian=False):
     omega_r = max( a_rj.abs() for a_rj in A1[:r,r].list() )
 
     if A1[0,0].abs()*omega_r >= alpha*(omega_1^2):
-        # Higham step (2)
-        # The top-left entry is our 1x1 pivot.
-        C = A1[1:n,0]
-        B = A1[1:,1:]
-
-        P2, L2, D2 = block_ldlt_naive(B - (C*C.transpose())/A1[0,0])
-
-        P1 = block_matrix(2,2, [[ZZ(1), ZZ(0)],
-                                [0*C,  P2]])
-        L1 = block_matrix(2,2, [[ZZ(1),                    ZZ(0)],
-                                [P2.transpose()*C/A1[0,0], L2]])
-        D1 = block_matrix(2,2, [[A1[0,0], ZZ(0)],
-                                [0*C,   D2]])
-
-        return (P1,L1,D1)
-
+        return pivot_one_by_one(A1)
 
     if A1[r,r].abs() > alpha*omega_r:
         # Higham step (3)
         # Another 1x1 pivot, but this time swapping indices 0,r.
-        P1 = copy(A1.matrix_space().identity_matrix())
-        P1.swap_rows(0,s)
-        A1 = P1.T * A1 * P1
-
-        # The top-left entry is now our 1x1 pivot.
-        C = A1[1:n,0]
-        B = A1[1:,1:]
-
-        P2, L2, D2 = block_ldlt_naive(B - (C*C.transpose())/A1[0,0])
-
-        P1 = P1*block_matrix(2,2, [[ZZ(1), ZZ(0)],
-                                   [0*C,   P2]])
-        L1 = block_matrix(2,2, [[ZZ(1),                    ZZ(0)],
-                                [P2.transpose()*C/A1[0,0], L2]])
-        D1 = block_matrix(2,2, [[A1[0,0], ZZ(0)],
-                                [0*C,     D2]])
-
-        return (P1,L1,D1)
+        return pivot_one_by_one(A1,r)
 
     # Higham step (4)
     # If we made it here, we have to do a 2x2 pivot.
-- 
2.44.2