# HG changeset patch
# User Jaroslav Hajek <highegg@gmail.com>
# Date 1259216522 -3600
# Node ID 87595d71400529ce0f39881563c5ad5d445e6437
# Parent  5b733adba096860741ac87f7fef9edbb8353e21e
move normest to linear-algebra, improve it

diff --git a/scripts/ChangeLog b/scripts/ChangeLog
--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -1,3 +1,8 @@
+2009-11-26  Jaroslav Hajek  <highegg@gmail.com>
+
+	* sparse/normest.m: Move to linear-algebra.
+	* linear-algebra/normest.m: Simplify. Don't form A'*A explicitly.
+
 2009-11-25  Jaroslav Hajek  <highegg@gmail.com>
 
 	* linear-algebra/isdefinite.m: Use Cholesky factorization.
diff --git a/scripts/linear-algebra/module.mk b/scripts/linear-algebra/module.mk
--- a/scripts/linear-algebra/module.mk
+++ b/scripts/linear-algebra/module.mk
@@ -15,6 +15,7 @@
   linear-algebra/krylov.m \
   linear-algebra/krylovb.m \
   linear-algebra/logm.m \
+  linear-algebra/normest.m \
   linear-algebra/null.m \
   linear-algebra/onenormest.m \
   linear-algebra/orth.m \
diff --git a/scripts/sparse/normest.m b/scripts/linear-algebra/normest.m
rename from scripts/sparse/normest.m
rename to scripts/linear-algebra/normest.m
--- a/scripts/sparse/normest.m
+++ b/scripts/linear-algebra/normest.m
@@ -1,4 +1,5 @@
 ## Copyright (C) 2006, 2007, 2008, 2009 David Bateman and Marco Caliari
+## Copyright (C) 2009 VZLU Prague
 ##
 ## This file is part of Octave.
 ##
@@ -28,60 +29,32 @@
 ## @code{normest} to converge.
 ## @end deftypefn
 
-function [e1, c] = normest (A, tol)
-  if (nargin < 2)
-    tol = 1e-6;
-  endif
-  if (isa (A, "single"))
-    if (tol < eps ("single"))
-      tol = eps ("single");
-    endif
-  else
-    if (tol < eps)
-      tol = eps
-    endif
-  endif
-  if (ndims (A) != 2)
+function [e, c] = normest (A, tol = 1e-6)
+  if (! ismatrix (A) || ndims (A) != 2)
     error ("normest: A must be a matrix");
   endif 
-  maxA = max (max (abs (A)));
+  if (! isfloat (A))
+    A = double (A);
+  endif
+  tol = max (tol, eps (class (A)));
   c = 0;
-  if (maxA == 0)
-    e1 = 0
-  else
+  x = norm (A, "columns").';
+  e = norm (x);
+  if (e > 0)
     [m, n] = size (A);
-    B = A / maxA;
-    Bt = B';
-    if (m > n)
-      tmp = B;
-      B = Bt;
-      Bt = tmp;
-    endif
+    x /= e;
     e0 = 0;
-    x = randn (min (m, n), 1);
-    e1 = norm (x);
-    x = x / e1;
-    e1 = sqrt (e1);
-    if (issparse (A))
-      while (abs (e1 - e0) > tol * e1)
-	e0 = e1;
-	x = B * (Bt * x);
-	e1 = norm (x);
-	x = x / e1;
-	e1 = sqrt (e1);
-	c = c + 1;
-      endwhile
-    else
-      B = B * Bt;
-      while (abs (e1 - e0) > tol * e1)
-	e0 = e1;
-	x = B * x;
-	e1 = norm (x);
-	x = x / e1;
-	e1 = sqrt (e1);
-	c = c + 1;
-      endwhile
-    endif
-    e1 = e1 * maxA;
+    while (abs (e - e0) > tol * e)
+      e0 = e;
+      y = A*x;
+      e = norm (y);
+      if (e == 0)
+        x = rand (n, 1);
+      else
+        x = A'*(y / e);
+      endif
+      x /= norm (x);
+      c += 1;
+    endwhile
   endif
 endfunction
diff --git a/scripts/sparse/module.mk b/scripts/sparse/module.mk
--- a/scripts/sparse/module.mk
+++ b/scripts/sparse/module.mk
@@ -7,7 +7,6 @@
   sparse/etreeplot.m \
   sparse/gplot.m \
   sparse/nonzeros.m \
-  sparse/normest.m \
   sparse/pcg.m \
   sparse/pcr.m \
   sparse/spalloc.m \