+ # The top-left 2x2 submatrix (starting at position k,k) is now
+ # our pivot.
+ E = A[k:k+2,k:k+2]
+ d.append(E)
+
+ C = A[k+2:n,k:k+2]
+ B = A[k+2:,k+2:]
+
+ # 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.
+ #
+ # The latter can be found much more 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())
+
+ # 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+j,k+2+i] = A[k+2+j,k+2+i] - schur_complement[j,i]
+ A[k+2+i,k+2+j] = A[k+2+j,k+2+i] # keep it symmetric!
+
+ # The on- and above-diagonal entries of "L" will be fixed
+ # later, so we only need to worry about the lower-left entry
+ # of the 2x2 identity matrix that belongs at the top of the
+ # new column of "L".
+ A[k+1,k] = 0
+ 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.
+ A[k+i+2,k+j] = CE_inverse[i,j]
+
+
+ k += 2