# HG changeset patch
# User Jaroslav Hajek <highegg@gmail.com>
# Date 1259216375 -3600
# Node ID 5b733adba096860741ac87f7fef9edbb8353e21e
# Parent  ecd750d1eabdc60a67c80ded35f1aa6ef955ff0e
base isdefinite on cholesky decomposition

diff --git a/scripts/ChangeLog b/scripts/ChangeLog
--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -1,3 +1,7 @@
+2009-11-25  Jaroslav Hajek  <highegg@gmail.com>
+
+	* linear-algebra/isdefinite.m: Use Cholesky factorization.
+
 2009-11-24  Jaroslav Hajek  <highegg@gmail.com>
 
 	* general/issymmetric.m: Move to linear-algebra.
diff --git a/scripts/linear-algebra/isdefinite.m b/scripts/linear-algebra/isdefinite.m
--- a/scripts/linear-algebra/isdefinite.m
+++ b/scripts/linear-algebra/isdefinite.m
@@ -21,7 +21,8 @@
 ## Return 1 if @var{x} is symmetric positive definite within the
 ## tolerance specified by @var{tol} or 0 if @var{x} is symmetric
 ## positive semidefinite.  Otherwise, return -1.  If @var{tol}
-## is omitted, use a tolerance equal to 100 times the machine precision.
+## is omitted, use a tolerance equal to 100 times the machine precision,
+## multiplied by the Frobeniusm norm of @var{x}.
 ## @seealso{issymmetric}
 ## @end deftypefn
 
@@ -31,30 +32,33 @@
 
 function retval = isdefinite (x, tol)
 
-  if (nargin == 1 || nargin == 2)
-    if (nargin == 1)
-      if (isa (x, "single"))
-	tol = 100 * eps("single");
-      else
-	tol = 100*eps; 
-      endif
-    endif
-    sym = ishermitian (x);
-    if (sym)
-      ## Matrix is symmetric, check eigenvalues.
-      mineig = min (eig (x));
-      if (mineig > tol)
-	retval = 1;
-      elseif (mineig > -tol)
-	retval = 0;
-      else
-	retval = -1;
-      endif
-    else
-      error ("isdefinite: matrix x must be symmetric");
-    endif
-  else
+  if (nargin < 1 || nargin > 2)
     print_usage ();
   endif
 
+  if (! ishermitian (x))
+    error ("isdefinite: x must be a hermitian matrix");
+  endif
+
+  if (! isfloat (x))
+    x = double (x);
+  endif
+
+  if (nargin == 1)
+    tol = 100 * eps(class (x)) * norm (x, "fro");
+  endif
+
+  e = tol * eye (rows (x));
+  [r, p] = chol (x - e);
+  if (p == 0)
+    retval = 1;
+  else
+    [r, p] = chol (x + e);
+    if (p == 0)
+      retval = 0;
+    else
+      retval = -1;
+    endif
+  endif
+
 endfunction