mjo/ldlt.py: make block_ldlt() work on Hermitian matrices, too.

diff --git a/mjo/ldlt.py b/mjo/ldlt.py
index 09b3bb008263281e9cc787764972bb14d40392e4..729e01152033772217910b1fc46adf50f5d54b8c 100644 (file)
--- a/mjo/ldlt.py
+++ b/mjo/ldlt.py
@@ -435,15 +435,13 @@ def block_ldlt(A):
         for i in range(n-k-1):
             for j in range(i+1):
                 A[k+1+i,k+1+j] = ( A[k+1+i,k+1+j] -
-                                   A[k,k+1+i]*A[k,k+1+j]/A[k,k] )
-                A[k+1+j,k+1+i] = A[k+1+i,k+1+j] # keep it symmetric!
+                                   A[k+1+i,k]*A[k,k+1+j]/A[k,k] )
+                A[k+1+j,k+1+i] = A[k+1+i,k+1+j].conjugate() # stay hermitian!
 
         for i in range(n-k-1):
             # Store the new (kth) column of "L" within the lower-
-            # left-hand corner of "A", being sure to set the lower-
-            # left entries from the upper-right ones to avoid
-            # collisions.
-            A[k+i+1,k] = A[k,k+1+i]/A[k,k]
+            # left-hand corner of "A".
+            A[k+i+1,k] /= A[k,k]
 
         # No return value, only the desired side effects of updating
         # p, d, and A.
@@ -460,7 +458,9 @@ def block_ldlt(A):
         if k == (n-1):
             # Handle this trivial case manually, since otherwise the
             # algorithm's references to the e.g. "subdiagonal" are
-            # meaningless.
+            # meaningless. The corresponding entry of "L" will be
+            # fixed later (since it's an on-diagonal element, it gets
+            # set to one eventually).
             d.append( matrix(ring, 1, [[A[k,k]]]) )
             k += 1
             continue
@@ -469,10 +469,8 @@ def block_ldlt(A):
         # kth column. This occurs prior to Step (1) in Higham,
         # but is part of Step (1) in Bunch and Kaufman. We adopt
         # Higham's "omega" notation instead of B&K's "lambda"
-        # because "lambda" can lead to some confusion. Beware:
-        # the subdiagonals of our matrix are being overwritten!
-        # So we actually use the corresponding row entries instead.
-        column_1_subdiag = [ a_ki.abs() for a_ki in A[k,k+1:].list() ]
+        # because "lambda" can lead to some confusion.
+        column_1_subdiag = [ a_ki.abs() for a_ki in A[k+1:,k].list() ]
         omega_1 = max([ a_ki for a_ki in column_1_subdiag ])
 
         if omega_1 == 0:
@@ -482,7 +480,9 @@ def block_ldlt(A):
             #   [ 0 B ]
             #
             # and we can simply skip to the next step after recording
-            # the 1x1 pivot "1" in the top-left position.
+            # the 1x1 pivot "a" in the top-left position. The entry "a"
+            # will be adjusted to "1" later on to ensure that "L" is
+            # (block) unit-lower-triangular.
             d.append( matrix(ring, 1, [[A[k,k]]]) )
             k += 1
             continue
@@ -508,9 +508,8 @@ def block_ldlt(A):
         # B&K's Step (3) where we find the largest off-diagonal entry
         # (in magniture) in column "r". Since the matrix is Hermitian,
         # we need only look at the above-diagonal entries to find the
-        # off-diagonal of maximal magnitude. (Beware: the subdiagonal
-        # entries are being overwritten.)
-        omega_r = max( a_rj.abs() for a_rj in A[:r,r].list() )
+        # off-diagonal of maximal magnitude.
+        omega_r = max( a_rj.abs() for a_rj in A[r,k:r].list() )
 
         if A[k,k].abs()*omega_r >= alpha*(omega_1**2):
             # Step (2) in Higham or Step (4) in B&K.
@@ -540,22 +539,22 @@ def block_ldlt(A):
         # We don't actually need the inverse of E, what we really need
         # is C*E.inverse(), and that can be found by setting
         #
-        #   C*E.inverse() == X   <====>   XE == C.
+        #   X = C*E.inverse()   <====>   XE = C.
         #
-        # The latter can be found much more easily by solving a system.
-        # Note: I do not actually know that sage solves the system more
+        # Then "X" can be found easily by solving a system.  Note: I
+        # do not actually know that sage solves the system more
         # intelligently, but this is still The Right Thing To Do.
         CE_inverse = E.solve_left(C)
 
-        schur_complement = B - (CE_inverse*C.transpose())
+        schur_complement = B - (CE_inverse*C.conjugate_transpose())
 
         # Compute the Schur complement that we'll work on during
         # the following iteration, and store it back in the lower-
         # right-hand corner of "A".
         for i in range(n-k-2):
             for j in range(i+1):
-                A[k+2+i,k+2+j] = A[k+2+i,k+2+j] - schur_complement[i,j]
-                A[k+2+j,k+2+i] = A[k+2+j,k+2+i] - schur_complement[j,i]
+                A[k+2+i,k+2+j] = schur_complement[i,j]
+                A[k+2+j,k+2+i] = schur_complement[j,i]
 
         # The on- and above-diagonal entries of "L" will be fixed
         # later, so we only need to worry about the lower-left entry
@@ -565,9 +564,7 @@ def block_ldlt(A):
         for i in range(n-k-2):
             for j in range(2):
                 # Store the new (k and (k+1)st) columns of "L" within
-                # the lower-left-hand corner of "A", being sure to set
-                # the lower-left entries from the upper-right ones to
-                # avoid collisions.
+                # the lower-left-hand corner of "A".
                 A[k+i+2,k+j] = CE_inverse[i,j]