- n = A.nrows()
-
- # Use the fraction field of the given matrix so that division will work
- # when (for example) our matrix consists of integer entries.
- ring = A.base_ring().fraction_field()
-
- if n == 0 or n == 1:
- # We can get n == 0 if someone feeds us a trivial matrix.
- P = matrix.identity(ring, n)
- L = matrix.identity(ring, n)
- D = A
- return (P,L,D)
-
- A1 = A.change_ring(ring)
- diags = A1.diagonal()
- s = diags.index(max(diags))
- P1 = copy(A1.matrix_space().identity_matrix())
- A1 = P1.T * A1 * P1
- alpha1 = A1[0,0]
-
- # Golub and Van Loan mention in passing what to do here. This is
- # only sensible if the matrix is positive-semidefinite, because we
- # are assuming that we can set everything else to zero as soon as
- # we hit the first on-diagonal zero.
- if alpha1 == 0:
- P = A1.matrix_space().identity_matrix()
- L = P
- D = A1.matrix_space().zero()
- return (P,L,D)
-
- v1 = A1[1:n,0]
- A2 = A1[1:,1:]
-
- P2, L2, D2 = ldlt_naive(A2 - (v1*v1.transpose())/alpha1)
-
- P1 = P1*block_matrix(2,2, [[ZZ(1), ZZ(0)],
- [0*v1, P2]])
- L1 = block_matrix(2,2, [[ZZ(1), ZZ(0)],
- [P2.transpose()*v1/alpha1, L2]])
- D1 = block_matrix(2,2, [[alpha1, ZZ(0)],
- [0*v1, D2]])
-
- return (P1,L1,D1)