mjo/ldlt.py: fix two bugs in the imperative block_ldlt() function.

diff --git a/mjo/ldlt.py b/mjo/ldlt.py
index 6b9908341add2dfaadb99c17a76dc7c4c19913dc..c8483c707afa6f3b67935e451b7c979401661b09 100644 (file)
--- a/mjo/ldlt.py
+++ b/mjo/ldlt.py
@@ -409,28 +409,28 @@ def block_ldlt(A):
             p[k] = p[s]
             p[s] = p_k
 
-            # Now the pivot is in the (k,k)th position.
-            d.append( matrix(ring, 1, [[A[k,k]]]) )
-
-            # 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-1):
-                for j in range(i+1):
-                    A[k+1+j,k+1+i] = ( A[k+1+j,k+1+i] -
-                                       A[k,k+1+j]*A[k,k+1+i]/alpha )
-                    A[k+1+i,k+1+j] = A[k+1+j,k+1+i] # keep it symmetric!
-
-            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]/alpha
-
-            # No return value, only the desired side effects of updating
-            # p, d, and A.
-            return
+        # Now the pivot is in the (k,k)th position.
+        d.append( matrix(ring, 1, [[A[k,k]]]) )
+
+        # 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-1):
+            for j in range(i+1):
+                A[k+1+j,k+1+i] = ( A[k+1+j,k+1+i] -
+                                   A[k,k+1+j]*A[k,k+1+i]/A[k,k] )
+                A[k+1+i,k+1+j] = A[k+1+j,k+1+i] # keep it symmetric!
+
+        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]
+
+        # No return value, only the desired side effects of updating
+        # p, d, and A.
+        return
 
     k = 0
     while k < n: