Mercurial > hg > octave-terminal
changeset 9871:87595d714005
move normest to linear-algebra, improve it
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Thu, 26 Nov 2009 07:22:02 +0100 |
| parents | 5b733adba096 |
| children | 72d6e0de76c7 |
| files | scripts/ChangeLog scripts/linear-algebra/module.mk scripts/linear-algebra/normest.m scripts/sparse/module.mk scripts/sparse/normest.m |
| diffstat | 4 files changed, 29 insertions(+), 51 deletions(-) [+] |
line wrap: on
line diff
--- 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.
--- 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 \
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