from sage.all import *

def is_positive_semidefinite_naive(A):
    r"""
    A naive positive-semidefinite check that tests the eigenvalues for
    nonnegativity.  We follow the sage convention that positive
    (semi)definite matrices must be symmetric or Hermitian.

    SETUP::

        sage: from mjo.ldlt import is_positive_semidefinite_naive

    TESTS:

    The trivial matrix is vaciously positive-semidefinite::

        sage: A = matrix(QQ, 0)
        sage: A
        []
        sage: is_positive_semidefinite_naive(A)
        True

    """
    if A.nrows() == 0:
        return True # vacuously
    return A.is_hermitian() and all( v >= 0 for v in A.eigenvalues() )

def ldlt_naive(A):
    r"""
    Perform a pivoted `LDL^{T}` factorization of the Hermitian
    positive-semidefinite matrix `A`.

    This is a naive, recursive implementation that is inefficient due
    to Python's lack of tail-call optimization. The pivot strategy is
    to choose the largest diagonal entry of the matrix at each step,
    and to permute it into the top-left position. Ultimately this
    results in a factorization `A = PLDL^{T}P^{T}`, where `P` is a
    permutation matrix, `L` is unit-lower-triangular, and `D` is
    diagonal decreasing from top-left to bottom-right.

    ALGORITHM:

    The algorithm is based on the discussion in Golub and Van Loan, but with
    some "typos" fixed.

    OUTPUT:

    A triple `(P,L,D)` such that `A = PLDL^{T}P^{T}` and where,

      * `P` is a permutaiton matrix
      * `L` is unit lower-triangular
      * `D` is a diagonal matrix whose entries are decreasing from top-left
        to bottom-right

    SETUP::

        sage: from mjo.ldlt import ldlt_naive, is_positive_semidefinite_naive

    EXAMPLES:

    All three factors should be the identity when the original matrix is::

        sage: I = matrix.identity(QQ,4)
        sage: P,L,D = ldlt_naive(I)
        sage: P == I and L == I and D == I
        True

    TESTS:

    Ensure that a "random" positive-semidefinite matrix is factored correctly::

        sage: set_random_seed()
        sage: n = ZZ.random_element(5)
        sage: A = matrix.random(QQ, n)
        sage: A = A*A.transpose()
        sage: is_positive_semidefinite_naive(A)
        True
        sage: P,L,D = ldlt_naive(A)
        sage: A == P*L*D*L.transpose()*P.transpose()
        True

    """
    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())
    P1.swap_rows(0,s)
    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)