- 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/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-
- # diagonal entries to find the off-diagonal of maximal magnitude.
- 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):
- 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.
- return pivot_one_by_one(A1,r)
-
- # Higham step (4)
- # If we made it here, we have to do a 2x2 pivot.
- P1 = copy(A1.matrix_space().identity_matrix())
- P1.swap_rows(1,r)
- A1 = P1.T * A1 * P1
-
- # The top-left 2x2 submatrix is now our pivot.
- E = A1[:2,:2]
- C = A1[2:n,0:2]
- B = A1[2:,2:]
-
- if B.nrows() == 0:
- # We have a two-by-two matrix that we can do nothing
- # useful with.
- P = matrix.identity(ring, n)
- L = matrix.identity(ring, n)
- D = A1
- return (P,L,D)
-
- P2, L2, D2 = block_ldlt_naive(B - (C*E.inverse()*C.transpose()))
-
- 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*E.inverse(), L2]])
- D1 = block_diagonal_matrix(E,D2)
-
- return (P1,L1,D1)
-
-
-def block_ldlt(A):
- r"""
- Perform a block-`LDL^{T}` factorization of the Hermitian
- matrix `A`.
-
- OUTPUT:
-
- A triple `(P,L,D)` such that `A = PLDL^{T}P^{T}` and where,
-
- * `P` is a permutation matrix
- * `L` is unit lower-triangular
- * `D` is a block-diagonal matrix whose blocks are of size
- one or two.
- """
-
- # 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".
- 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 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.
-
- Note that ``A`` is passed in by reference here, so it doesn't
- matter if we shadow the name ``A`` with itself.
- """
- if s > k:
- # s == k would swap row/column k with itself, and we don't
- # actually want to perform the identity permutation.
- # We don't have to permute "L" separately so long as "L"
- # is stored within "A".
- A.swap_columns(k,s)
- A.swap_rows(k,s)
-
- # Update the permutation "matrix" with the swap we just did.
- p_k = p[k]
- 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]/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:
- # 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.
- 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. 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,1:].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 "1" in the top-left position.
- 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. (Beware: the subdiagonal
- # entries are being overwritten.)
- omega_r = max( a_rj.abs() for a_rj in A[:r,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(A1,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.
- 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)