X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Fldlt.py;h=1c2663ea3d574360059b07d22f9af8163ebdc2e7;hb=fe85934cbfb26ff0dc7e5ef81411fe73eed739ff;hp=a86dbbefe5f396a197a76af33dc60583cb7cb971;hpb=b1cae89cc2e8049c43638d7a8bed2256a72f1650;p=sage.d.git

diff --git a/mjo/ldlt.py b/mjo/ldlt.py
index a86dbbe..1c2663e 100644
--- a/mjo/ldlt.py
+++ b/mjo/ldlt.py
@@ -24,3 +24,50 @@ def is_positive_semidefinite_naive(A):
     if A.nrows() == 0:
         return True # vacuously
     return A.is_hermitian() and all( v >= 0 for v in A.eigenvalues() )
+
+
+def is_positive_semidefinite(A):
+    r"""
+    A fast positive-semidefinite check based on the block-LDLT
+    factorization.
+
+    SETUP::
+
+        sage: from mjo.ldlt import (is_positive_semidefinite,
+        ....:                       is_positive_semidefinite_naive)
+
+    TESTS:
+
+    Check that the naive and fast answers are the same, in general::
+
+        sage: set_random_seed()
+        sage: F = NumberField(x^2 + 1, 'I')
+        sage: from sage.misc.prandom import choice
+        sage: ring = choice([ZZ,QQ,F])
+        sage: A = matrix.random(ring, 10)
+        sage: is_positive_semidefinite(A) == is_positive_semidefinite_naive(A)
+        True
+
+    Check that the naive and fast answers are the same for a Hermitian
+    matrix::
+
+        sage: set_random_seed()
+        sage: F = NumberField(x^2 + 1, 'I')
+        sage: from sage.misc.prandom import choice
+        sage: ring = choice([ZZ,QQ,F])
+        sage: A = matrix.random(ring, 10); A = A + A.conjugate_transpose()
+        sage: is_positive_semidefinite(A) == is_positive_semidefinite_naive(A)
+        True
+
+    """
+    if not A.is_hermitian():
+        return False
+    _,_,d = A._block_ldlt()
+    for d_i in d:
+        if d_i.nrows() == 1:
+            if d_i < 0:
+                return False
+        else:
+            # A 2x2 block indicates that it's indefinite
+            return False
+    return True