Mercurial > hg > octave-lyh
annotate scripts/sparse/normest.m @ 8920:eb63fbe60fab
update copyright notices
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Sat, 07 Mar 2009 10:41:27 -0500 |
| parents | cadc73247d65 |
| children | 1bf0ce0930be |
| rev | line source |
|---|---|
| 8920 | 1 ## Copyright (C) 2006, 2007, 2008, 2009 David Bateman and Marco Caliari |
| 6212 | 2 ## |
| 7016 | 3 ## This file is part of Octave. |
| 6212 | 4 ## |
| 7016 | 5 ## Octave is free software; you can redistribute it and/or modify it |
| 6 ## under the terms of the GNU General Public License as published by | |
| 7 ## the Free Software Foundation; either version 3 of the License, or (at | |
| 8 ## your option) any later version. | |
| 9 ## | |
| 10 ## Octave is distributed in the hope that it will be useful, but | |
| 11 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
| 12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
| 13 ## General Public License for more details. | |
| 6212 | 14 ## |
| 15 ## You should have received a copy of the GNU General Public License | |
| 7016 | 16 ## along with Octave; see the file COPYING. If not, see |
| 17 ## <http://www.gnu.org/licenses/>. | |
| 6212 | 18 |
| 19 ## -*- texinfo -*- | |
| 20 ## @deftypefn {Function File} {[@var{n}, @var{c}] =} normest (@var{a}, @var{tol}) | |
| 6222 | 21 ## Estimate the 2-norm of the matrix @var{a} using a power series |
| 6212 | 22 ## analysis. This is typically used for large matrices, where the cost |
| 6222 | 23 ## of calculating the @code{norm (@var{a})} is prohibitive and an approximation |
| 6212 | 24 ## to the 2-norm is acceptable. |
| 25 ## | |
| 26 ## @var{tol} is the tolerance to which the 2-norm is calculated. By default | |
| 27 ## @var{tol} is 1e-6. @var{c} returns the number of iterations needed for | |
| 28 ## @code{normest} to converge. | |
| 29 ## @end deftypefn | |
| 30 | |
| 6498 | 31 function [e1, c] = normest (A, tol) |
| 6212 | 32 if (nargin < 2) |
| 33 tol = 1e-6; | |
| 34 endif | |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
35 if (isa (A, "single")) |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
36 if (tol < eps ("single")) |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
37 tol = eps ("single"); |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
38 endif |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
39 else |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
40 if (tol < eps) |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
41 tol = eps |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
42 endif |
| 6212 | 43 endif |
| 8507 | 44 if (ndims (A) != 2) |
| 6498 | 45 error ("normest: A must be a matrix"); |
| 6212 | 46 endif |
| 6498 | 47 maxA = max (max (abs (A))); |
| 6212 | 48 c = 0; |
| 49 if (maxA == 0) | |
| 50 e1 = 0 | |
| 51 else | |
| 6498 | 52 [m, n] = size (A); |
| 6212 | 53 B = A / maxA; |
| 54 Bt = B'; | |
| 55 if (m > n) | |
| 56 tmp = B; | |
| 57 B = Bt; | |
| 58 Bt = tmp; | |
| 59 endif | |
| 60 e0 = 0; | |
| 6498 | 61 x = randn (min (m, n), 1); |
| 62 e1 = norm (x); | |
| 6212 | 63 x = x / e1; |
| 6498 | 64 e1 = sqrt (e1); |
| 65 if (issparse (A)) | |
| 66 while (abs (e1 - e0) > tol * e1) | |
| 6212 | 67 e0 = e1; |
| 68 x = B * (Bt * x); | |
| 6498 | 69 e1 = norm (x); |
| 6212 | 70 x = x / e1; |
| 6498 | 71 e1 = sqrt (e1); |
| 6212 | 72 c = c + 1; |
| 73 endwhile | |
| 74 else | |
| 6213 | 75 B = B * Bt; |
| 6498 | 76 while (abs (e1 - e0) > tol * e1) |
| 6212 | 77 e0 = e1; |
| 78 x = B * x; | |
| 6498 | 79 e1 = norm (x); |
| 6212 | 80 x = x / e1; |
| 6498 | 81 e1 = sqrt (e1); |
| 6212 | 82 c = c + 1; |
| 83 endwhile | |
| 84 endif | |
| 85 e1 = e1 * maxA; | |
| 86 endif | |
| 87 endfunction |