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
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
|
function lintest (A,b)
%LINTEST test A*x=b, using linfactor, x=A\b, and (ack!) the explicit inv(A).
% The results printed include the breakeven point, which is the number of
% systems Ax=b that must be solved with the same A but different b for the
% inv(A) method to be faster than linfactor. Using inv(A) is always
% numerically dubious, and typically slower. However because of MATLAB's
% interpretive overhead, linfactor can be slightly slower in its
% forward/backsolves. A true linfactor mexFunction would probably be just as
% fast as inv(A)*b when A is full. You should never ever use inv(A) to solve
% a linear system.
%
% Example:
% load west0479
% b = rand (size (west0479,1), 1) ;
% lintest (west0479, b) ;
%
% See also linfactor, lintests, mldivide
% Copyright 2007, Timothy A. Davis
%-------------------------------------------------------------------------------
% warmup, to make sure functions are loaded, for accurate timings
%-------------------------------------------------------------------------------
F = linfactor (1) ;
x = linfactor (F, 1) ; %#ok
F = linfactor (sparse (1)) ;
x = linfactor (F, 1) ; %#ok
F = linfactor (sparse (-1)) ;
x = linfactor (F, 1) ; %#ok
S = inv (1) ;
x = S*1 ; %#ok
S = inv (sparse (1)) ;
x = S*1 ; %#ok
S = inv (sparse (-1)) ;
x = S*1 ; %#ok
x = rand (2) \ rand (2,1) ; %#ok
x = sparse (rand (2)) \ rand (2,1) ; %#ok
clear x F S
%-------------------------------------------------------------------------------
% linfactor
%-------------------------------------------------------------------------------
% do this several times for accurate timings
t1 = 0 ;
trials = 0 ;
while (t1 < 1)
[F, t] = linfactor (A) ; % factorize A
t1 = t1 + t ;
trials = trials + 1 ;
end
t1 = t1 / trials ;
% do this several times for accurate timings
t2 = 0 ;
trials = 0 ;
while (t2 < 1)
[x, t] = linfactor (F, b) ; % use the factors to solve Ax=b
t2 = t2 + t ;
trials = trials + 1 ;
end
t2 = t2 / trials ;
resid = norm (A*x-b,1) / (norm (A,1) * norm (x,1) + norm (b,1)) ;
fprintf ('%-16s factor time: %10.6f solve time: %10.6f resid: %8.2e\n', ...
F.kind (1:(find(F.kind == ':', 1, 'first'))), t1, t2, resid) ;
%-------------------------------------------------------------------------------
% inv
%-------------------------------------------------------------------------------
% Try again using inv(A)*b. This is a really horrible way to solve Ax=b. I'm
% doing it here precisely to show that it is typically slower, except when there
% are a huge number of right-hand-sides to solve and A is small. inv(A) also
% tends to be less accurate, but random matrices do not trigger that problem.
% It fails hopelessly when A is large, sparse, and where max(diff(r)) is large
% where [p,q,r,s] = dmperm(A). Never, ever use inv(A) to solve a linear system.
% Oh, did I tell you never to use inv(A) to solve Ax=b?
try
% do this several times for accurate timings
t3 = 0 ;
trials = 0 ;
while (t3 < 1)
tic ;
S = inv (A) ; %#ok
t = toc ;
t3 = t3 + t ;
trials = trials + 1 ;
end
t3 = t3 / trials ;
% do this several times for accurate timings
t4 = 0 ;
trials = 0 ;
while (t4 < 1)
tic ;
x = S*b ;
t = toc ;
t4 = t4 + t ;
trials = trials + 1 ;
end
t4 = t4 / trials ;
resid = norm (A*x-b,1) / (norm (A,1) * norm (x,1) + norm (b,1)) ;
fprintf ('%-16s factor time: %10.6f solve time: %10.6f resid: %8.2e\n', ...
'inv(A)', t3, t4, resid) ;
% determine the breakeven point where using inv(A) is faster
nrhs = max (1, ceil ((t3 - t1) / (t2 - t4))) ;
if (t1 < t3 & t2 < t4) %#ok
fprintf ('inv(A) breakeven: never\n') ;
elseif (t3 < t1 & t4 < t2) %#ok
fprintf ('inv(A) breakeven: > 0\n') ;
elseif (t1 < t3 & t4 < t2) %#ok
fprintf ('inv(A) breakeven: > %d\n', nrhs) ;
else
fprintf ('inv(A) breakeven: < %d\n', nrhs) ;
end
catch
% inv(A) probably ran out of memory
fprintf ('inv(A) failed\n') ;
fprintf ('inv(A) breakeven: never\n') ;
end
%-------------------------------------------------------------------------------
% backslash
%-------------------------------------------------------------------------------
% do this several times for accurate timings
t0 = 0 ;
trials = 0 ;
while (t0 < 1)
tic ;
x = A\b ; % solve Ax=b
t = toc ;
t0 = t0 + t ;
trials = trials + 1 ;
end
t0 = t0 / trials ;
resid = norm (A*x-b,1) / (norm (A,1) * norm (x,1) + norm (b,1)) ;
fprintf ('%-16s total time: %10.6f resid: %8.2e\n', ...
'x = A\b', t0, resid) ;
fprintf ('\n') ;
|
