--- a/scripts/ChangeLog	Thu Nov 26 07:19:35 2009 +0100
+++ b/scripts/ChangeLog	Thu Nov 26 07:22:02 2009 +0100
@@ -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.

--- a/scripts/linear-algebra/module.mk	Thu Nov 26 07:19:35 2009 +0100
+++ b/scripts/linear-algebra/module.mk	Thu Nov 26 07:22:02 2009 +0100
@@ -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 \

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/linear-algebra/normest.m	Thu Nov 26 07:22:02 2009 +0100
@@ -0,0 +1,60 @@
+## Copyright (C) 2006, 2007, 2008, 2009 David Bateman and Marco Caliari
+## Copyright (C) 2009 VZLU Prague
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{n}, @var{c}] =} normest (@var{a}, @var{tol})
+## Estimate the 2-norm of the matrix @var{a} using a power series
+## analysis.  This is typically used for large matrices, where the cost
+## of calculating the @code{norm (@var{a})} is prohibitive and an approximation
+## to the 2-norm is acceptable.
+##
+## @var{tol} is the tolerance to which the 2-norm is calculated.  By default
+## @var{tol} is 1e-6.  @var{c} returns the number of iterations needed for
+## @code{normest} to converge.
+## @end deftypefn
+
+function [e, c] = normest (A, tol = 1e-6)
+  if (! ismatrix (A) || ndims (A) != 2)
+    error ("normest: A must be a matrix");
+  endif 
+  if (! isfloat (A))
+    A = double (A);
+  endif
+  tol = max (tol, eps (class (A)));
+  c = 0;
+  x = norm (A, "columns").';
+  e = norm (x);
+  if (e > 0)
+    [m, n] = size (A);
+    x /= e;
+    e0 = 0;
+    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

--- a/scripts/sparse/module.mk	Thu Nov 26 07:19:35 2009 +0100
+++ b/scripts/sparse/module.mk	Thu Nov 26 07:22:02 2009 +0100
@@ -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 \

--- a/scripts/sparse/normest.m	Thu Nov 26 07:19:35 2009 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,87 +0,0 @@
-## Copyright (C) 2006, 2007, 2008, 2009 David Bateman and Marco Caliari
-##
-## This file is part of Octave.
-##
-## Octave is free software; you can redistribute it and/or modify it
-## under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 3 of the License, or (at
-## your option) any later version.
-##
-## Octave is distributed in the hope that it will be useful, but
-## WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-## General Public License for more details.
-##
-## You should have received a copy of the GNU General Public License
-## along with Octave; see the file COPYING.  If not, see
-## <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn {Function File} {[@var{n}, @var{c}] =} normest (@var{a}, @var{tol})
-## Estimate the 2-norm of the matrix @var{a} using a power series
-## analysis.  This is typically used for large matrices, where the cost
-## of calculating the @code{norm (@var{a})} is prohibitive and an approximation
-## to the 2-norm is acceptable.
-##
-## @var{tol} is the tolerance to which the 2-norm is calculated.  By default
-## @var{tol} is 1e-6.  @var{c} returns the number of iterations needed for
-## @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)
-    error ("normest: A must be a matrix");
-  endif 
-  maxA = max (max (abs (A)));
-  c = 0;
-  if (maxA == 0)
-    e1 = 0
-  else
-    [m, n] = size (A);
-    B = A / maxA;
-    Bt = B';
-    if (m > n)
-      tmp = B;
-      B = Bt;
-      Bt = tmp;
-    endif
-    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;
-  endif
-endfunction