Mercurial > hg > octave-lojdl
changeset 6212:cc34c0be7b00
[project @ 2006-12-08 11:17:12 by dbateman]
author | dbateman |
---|---|
date | Fri, 08 Dec 2006 11:18:15 +0000 |
parents | 778fd090e505 |
children | 0a259ae4375e |
files | scripts/ChangeLog scripts/sparse/normest.m |
diffstat | 2 files changed, 86 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,7 @@ +2006-12-08 David Bateman <dbateman@free.fr> + + * sparse/normest.m: New file. + 2006-12-06 Michael Goffioul <michael.goffioul@swing.be>. * miscellaneous/copyfile.m, miscellaneous/movefile.m:
new file mode 100644 --- /dev/null +++ b/scripts/sparse/normest.m @@ -0,0 +1,82 @@ +## 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} {[@var{n}, @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 @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 (tol < eps) + tol = eps + endif + if (ndims(A) != 2) + error("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