Mercurial > hg > octave-nkf

comparison scripts/sparse/normest.m @ 6498:2c85044aa63f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-04-05 17:59:47 by jwe]
author jwe
date Thu, 05 Apr 2007 17:59:47 +0000
parents 07d967f75dba
children 93c65f2a5668
comparison
equal deleted inserted replaced
6497:fc8ed0c77e08 6498:2c85044aa63f
25 ## @var{tol} is the tolerance to which the 2-norm is calculated. By default 25 ## @var{tol} is the tolerance to which the 2-norm is calculated. By default
26 ## @var{tol} is 1e-6. @var{c} returns the number of iterations needed for 26 ## @var{tol} is 1e-6. @var{c} returns the number of iterations needed for
27 ## @code{normest} to converge. 27 ## @code{normest} to converge.
28 ## @end deftypefn 28 ## @end deftypefn
29 29
30 function [e1, c] = normest(A, tol) 30 function [e1, c] = normest (A, tol)
31 if (nargin < 2) 31 if (nargin < 2)
32 tol = 1e-6; 32 tol = 1e-6;
33 endif 33 endif
34 if (tol < eps) 34 if (tol < eps)
35 tol = eps 35 tol = eps
36 endif 36 endif
37 if (ndims(A) != 2) 37 if (ndims(A) != 2)
38 error("A must be a matrix"); 38 error ("normest: A must be a matrix");
39 endif 39 endif
40 maxA = max(max(abs(A))); 40 maxA = max (max (abs (A)));
41 c = 0; 41 c = 0;
42 if (maxA == 0) 42 if (maxA == 0)
43 e1 = 0 43 e1 = 0
44 else 44 else
45 [m, n] = size(A); 45 [m, n] = size (A);
46 B = A / maxA; 46 B = A / maxA;
47 Bt = B'; 47 Bt = B';
48 if (m > n) 48 if (m > n)
49 tmp = B; 49 tmp = B;
50 B = Bt; 50 B = Bt;
51 Bt = tmp; 51 Bt = tmp;
52 endif 52 endif
53 e0 = 0; 53 e0 = 0;
54 x = randn(min(m,n),1); 54 x = randn (min (m, n), 1);
55 e1 = norm(x); 55 e1 = norm (x);
56 x = x / e1; 56 x = x / e1;
57 e1 = sqrt(e1); 57 e1 = sqrt (e1);
58 if (issparse(A)) 58 if (issparse (A))
59 while (abs(e1 - e0) > tol * e1) 59 while (abs (e1 - e0) > tol * e1)
60 e0 = e1; 60 e0 = e1;
61 x = B * (Bt * x); 61 x = B * (Bt * x);
62 e1 = norm(x); 62 e1 = norm (x);
63 x = x / e1; 63 x = x / e1;
64 e1 = sqrt(e1); 64 e1 = sqrt (e1);
65 c = c + 1; 65 c = c + 1;
66 endwhile 66 endwhile
67 else 67 else
68 B = B * Bt; 68 B = B * Bt;
69 while (abs(e1 - e0) > tol * e1) 69 while (abs (e1 - e0) > tol * e1)
70 e0 = e1; 70 e0 = e1;
71 x = B * x; 71 x = B * x;
72 e1 = norm(x); 72 e1 = norm (x);
73 x = x / e1; 73 x = x / e1;
74 e1 = sqrt(e1); 74 e1 = sqrt (e1);
75 c = c + 1; 75 c = c + 1;
76 endwhile 76 endwhile
77 endif 77 endif
78 e1 = e1 * maxA; 78 e1 = e1 * maxA;
79 endif 79 endif