1function est = mynormest1 (L, U, P, Q)
2%MYNORMEST1 estimate norm(A,1), using LU factorization (L*U = P*A*Q).
3%
4% Example:
5%   est = mynormest1 (L, U, P, Q)
6% See also: testall
7
8% Copyright 2006-2012, Timothy A. Davis, http://www.suitesparse.com
9
10n = size (L,1) ;
11est = 0 ;
12S = zeros (n,1) ;
13
14for k = 1:5
15
16    if k == 1
17        x = ones (n,1) / n ;
18    else
19
20        j = find (abs (x) == max (abs (x))) ;
21        j = j (1) ;
22        x = zeros (n,1) ;
23        x (j) = 1 ;
24
25        % fprintf ('eka: k %d j %d est %g\n', k, j, est) ;
26    end
27
28
29    % x=A\x, but use the existing P*A*Q=L*U factorization
30
31    x = Q * (U \ (L \ (P*x))) ;
32
33    est_old = est ;
34    est = sum (abs (x)) ;
35
36    unchanged = 1 ;
37    for i = 1:n
38        if (x (i) >= 0)
39            s = 1 ;
40        else
41            s = -1 ;
42        end
43        if (s ~= S (i))
44            S (i) = s ;
45            unchanged = 0 ;
46        end
47    end
48
49    if (any (S ~= signum (x)))
50        S'                                                                  %#ok
51        signum(x)'                                                          %#ok
52        error ('Hey!') ;
53    end
54
55    if k > 1 & (est <= est_old | unchanged)                                 %#ok
56        break ;
57    end
58    x = S ;
59
60    % x=A'\x, but use the existing P*A*Q=L*U factorization
61    x = P' * (L' \ (U' \ (Q'*x))) ;
62
63    if k > 1
64        jnew = find (abs (x) == max (abs (x))) ;
65        if (jnew == j)
66            break ;
67        end
68    end
69
70end
71
72for k = 1:n
73    x (k) = power (-1, k+1) * (1 + ((k-1)/(n-1))) ;
74end
75
76% x=A\x, but use the existing P*A*Q=L*U factorization
77x = Q * (U \ (L \ (P*x))) ;
78
79est_new = 2 * sum (abs (x)) / (3 * n) ;
80if (est_new > est)
81    est = est_new ;
82end
83
84
85
86
87function s = signum (x)
88%SIGNUM compute sign of x
89s = ones (length (x),1) ;
90s (find (x < 0)) = -1 ;     %#ok
91