Mercurial > octave-nkf

changeset 9870:5b733adba096

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
base isdefinite on cholesky decomposition
author Jaroslav Hajek <highegg@gmail.com>
date Thu, 26 Nov 2009 07:19:35 +0100
parents ecd750d1eabd
children 87595d714005
files scripts/ChangeLog scripts/linear-algebra/isdefinite.m
diffstat 2 files changed, 32 insertions(+), 24 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/ChangeLog	Tue Nov 24 15:02:04 2009 +0100
+++ b/scripts/ChangeLog	Thu Nov 26 07:19:35 2009 +0100
@@ -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.
--- a/scripts/linear-algebra/isdefinite.m	Tue Nov 24 15:02:04 2009 +0100
+++ b/scripts/linear-algebra/isdefinite.m	Thu Nov 26 07:19:35 2009 +0100
@@ -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