New version of normest.

[David Bateman]

Ok, thanks to an off-line discussion with Marco Caliari, he helped me
rewrite the normest code to simplify the kernel of the loop, fix the
exit condition and use a random starting vector as something that is
less likely to get stuck with incorrect norm estimates. His example was
normest(toeplitz([-2 1 0 0],[-2 1 0 0])) where matlab, and my previous
code gave the second largest singular value as the 2-norm estimate due
to the starting vector being sum(A*A')'. He also suggested some other
improvements, including Aitken's acceleration for the power series, but
in the benchmarking I tried it didn't seem to offer any improvements. In
any case the attached code now takes the same number of iterations as
matlab's code offers better convergence and takes significantly less
time if matlab's normest and this code are both run under octave. I
supply the code as a patch.

Regards
David

*** ./scripts/sparse/normest.m.orig2    2006-12-06 18:11:40.100365205 +0100
--- ./scripts/sparse/normest.m  2006-12-06 21:53:03.057765160 +0100
***************
*** 0 ****
--- 1,57 ----
+ ## Copyright (C) 2006 David Bateman
+ ## Copyright (C) 2006 Marco Caliari
+ ##
+ ## This program 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 2 of the License, or
+ ## (at your option) any later version.
+ ##
+ ## This program 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 this program; if not, write to the Free Software
+ ## Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+ ## 02110-1301  USA
+ 
+ ## -*- texinfo -*-
+ ## @deftypefn {Function File} address@hidden, @var{c}] =} normest (@var{a}, 
@var{tol})
+ ##
+ ## Estimates 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 @address@hidden)} 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;
+   end
+   if (tol < eps)
+     tol = eps
+   end
+   if (ndims(A) ~= 2)
+     error('A must be a matrix');
+   end
+   B = A*A';
+   c = 0;
+   x = randn(size(A,2),1);
+   e0 = 0;
+   e1 = norm(x);
+   x = x / e1;
+   e1 = sqrt(e1);
+   while (abs(e1 - e0) > tol * e1)
+     e0 = e1;
+     x = B * x;
+     e1 = norm(x);
+     x = x / e1;
+     e1 = sqrt(e1);
+     c = c+1;
+   end
+ end