-
- # We have to make at least one copy of the input matrix so that we
- # can change the base ring to its fraction field. Both "L" and the
- # intermediate Schur complements will potentially have entries in
- # the fraction field. However, we don't need to make *two* copies.
- # We can't store the entries of "D" and "L" in the same matrix if
- # "D" will contain any 2x2 blocks; but we can still store the
- # entries of "L" in the copy of "A" that we're going to make.
- # Contrast this with the non-block LDL^T factorization where the
- # entries of both "L" and "D" overwrite the lower-left half of "A".
- #
- # This grants us an additional speedup, since we don't have to
- # permute the rows/columns of "L" *and* "A" at each iteration.
- ring = A.base_ring().fraction_field()
- A = A.change_ring(ring)
- MS = A.matrix_space()
-
- # The magic constant used by Bunch-Kaufman
- alpha = (1 + ZZ(17).sqrt()) * ~ZZ(8)
-
- # Keep track of the permutations and diagonal blocks in a vector
- # rather than in a matrix, for efficiency.
- n = A.nrows()
- p = list(range(n))
- d = []
-
- def swap_rows_columns(M, k, s):
- r"""
- Swap rows/columns ``k`` and ``s`` of the matrix ``M``, and update
- the list ``p`` accordingly.
- """
- if s > k:
- # s == k would swap row/column k with itself, and we don't
- # actually want to perform the identity permutation. If
- # you work out the recursive factorization by hand, you'll
- # notice that the rows/columns of "L" need to be permuted
- # as well. A nice side effect of storing "L" within "A"
- # itself is that we can skip that step. The first column
- # of "L" is hit by all of the transpositions in
- # succession, and the second column is hit by all but the
- # first transposition, and so on.
- M.swap_columns(k,s)
- M.swap_rows(k,s)
-
- p_k = p[k]
- p[k] = p[s]
- p[s] = p_k
-
- # No return value, we're only interested in the "side effects"
- # of modifing the matrix M (by reference) and the permutation
- # list p (which is in scope when this function is defined).
- return
-
-
- def pivot1x1(M, k, s):
- r"""
- Perform a 1x1 pivot swapping rows/columns `k` and `s >= k`.
- Relies on the fact that matrices are passed by reference,
- since for performance reasons this routine should overwrite
- its argument. Updates the local variables ``p`` and ``d`` as
- well.
- """
- swap_rows_columns(M,k,s)
-
- # 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+i,k+1+j] = ( A[k+1+i,k+1+j] -
- 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".
- A[k+i+1,k] /= A[k,k]
-
- # No return value, only the desired side effects of updating
- # p, d, and A.
- return
-
- k = 0
- while k < n:
- # At each step, we're considering the k-by-k submatrix
- # contained in the lower-right half of "A", because that's
- # where we're storing the next iterate. So our indices are
- # always "k" greater than those of Higham or B&K. Note that
- # ``n == 0`` is handled by skipping this loop entirely.
-
- if k == (n-1):
- # Handle this trivial case manually, since otherwise the
- # algorithm's references to the e.g. "subdiagonal" are
- # 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
-
- # Find the largest subdiagonal entry (in magnitude) in the
- # 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.
- 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:
- # In this case, our matrix looks like
- #
- # [ a 0 ]
- # [ 0 B ]
- #
- # and we can simply skip to the next step after recording
- # 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
-
- if A[k,k].abs() > alpha*omega_1:
- # This is the first case in Higham's Step (1), and B&K's
- # Step (2). Note that we have skipped the part of B&K's
- # Step (1) where we determine "r", since "r" is not yet
- # needed and we may waste some time computing it
- # otherwise. We are performing a 1x1 pivot, but the
- # rows/columns are already where we want them, so nothing
- # needs to be permuted.
- pivot1x1(A,k,k)
- k += 1
- continue
-
- # Now back to Step (1) of Higham, where we find the index "r"
- # that corresponds to omega_1. This is the "else" branch of
- # Higham's Step (1).
- r = k + 1 + column_1_subdiag.index(omega_1)
-
- # Continuing the "else" branch of Higham's Step (1), and onto
- # 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.
- 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.
- pivot1x1(A,k,k)
- k += 1
- continue
-
- if A[r,r].abs() > alpha*omega_r:
- # This is Step (3) in Higham or Step (5) in B&K. Still a 1x1
- # pivot, but this time we need to swap rows/columns k and r.
- pivot1x1(A,k,r)
- k += 1
- continue
-
- # If we've made it this far, we're at Step (4) in Higham or
- # Step (6) in B&K, where we perform a 2x2 pivot.
- swap_rows_columns(A,k+1,r)
-
- # 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
- #
- # X = C*E.inverse() <====> XE = C.
- #
- # 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.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] = 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
- # 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".
- A[k+i+2,k+j] = CE_inverse[i,j]
-
-
- k += 2
-
- MS = A.matrix_space()
- P = MS.matrix(lambda i,j: p[j] == i)
-
- # Warning: when n == 0, this works, but returns a matrix
- # whose (nonexistent) entries are in ZZ rather than in
- # the base ring of P and L.
- D = block_diagonal_matrix(d)
-
- # Overwrite the diagonal and upper-right half of "A",
- # since we're about to return it as the unit-lower-
- # triangular "L".
- for i in range(n):
- A[i,i] = 1
- for j in range(i+1,n):
- A[i,j] = 0
-
- return (P,A,D)