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