Mercurial > octave-nkf
view scripts/sparse/normest.m @ 6212:cc34c0be7b00
[project @ 2006-12-08 11:17:12 by dbateman]
author | dbateman |
---|---|
date | Fri, 08 Dec 2006 11:18:15 +0000 |
parents | |
children | 0a259ae4375e |
line wrap: on
line source
## 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