2
%LINTEST test A*x=b, using linfactor, x=A\b, and (ack!) the explicit inv(A).
3
% The results printed include the breakeven point, which is the number of
4
% systems Ax=b that must be solved with the same A but different b for the
5
% inv(A) method to be faster than linfactor. Using inv(A) is always
6
% numerically dubious, and typically slower. However because of MATLAB's
7
% interpretive overhead, linfactor can be slightly slower in its
8
% forward/backsolves. A true linfactor mexFunction would probably be just as
9
% fast as inv(A)*b when A is full. You should never ever use inv(A) to solve
14
% b = rand (size (west0479,1), 1) ;
15
% lintest (west0479, b) ;
17
% See also linfactor, lintests, mldivide
19
% Copyright 2007, Timothy A. Davis
21
%-------------------------------------------------------------------------------
22
% warmup, to make sure functions are loaded, for accurate timings
23
%-------------------------------------------------------------------------------
26
x = linfactor (F, 1) ; %#ok
27
F = linfactor (sparse (1)) ;
28
x = linfactor (F, 1) ; %#ok
29
F = linfactor (sparse (-1)) ;
30
x = linfactor (F, 1) ; %#ok
33
S = inv (sparse (1)) ;
35
S = inv (sparse (-1)) ;
37
x = rand (2) \ rand (2,1) ; %#ok
38
x = sparse (rand (2)) \ rand (2,1) ; %#ok
41
%-------------------------------------------------------------------------------
43
%-------------------------------------------------------------------------------
45
% do this several times for accurate timings
49
[F, t] = linfactor (A) ; % factorize A
55
% do this several times for accurate timings
59
[x, t] = linfactor (F, b) ; % use the factors to solve Ax=b
65
resid = norm (A*x-b,1) / (norm (A,1) * norm (x,1) + norm (b,1)) ;
67
fprintf ('%-16s factor time: %10.6f solve time: %10.6f resid: %8.2e\n', ...
68
F.kind (1:(find(F.kind == ':', 1, 'first'))), t1, t2, resid) ;
70
%-------------------------------------------------------------------------------
72
%-------------------------------------------------------------------------------
74
% Try again using inv(A)*b. This is a really horrible way to solve Ax=b. I'm
75
% doing it here precisely to show that it is typically slower, except when there
76
% are a huge number of right-hand-sides to solve and A is small. inv(A) also
77
% tends to be less accurate, but random matrices do not trigger that problem.
78
% It fails hopelessly when A is large, sparse, and where max(diff(r)) is large
79
% where [p,q,r,s] = dmperm(A). Never, ever use inv(A) to solve a linear system.
80
% Oh, did I tell you never to use inv(A) to solve Ax=b?
84
% do this several times for accurate timings
96
% do this several times for accurate timings
104
trials = trials + 1 ;
108
resid = norm (A*x-b,1) / (norm (A,1) * norm (x,1) + norm (b,1)) ;
110
fprintf ('%-16s factor time: %10.6f solve time: %10.6f resid: %8.2e\n', ...
111
'inv(A)', t3, t4, resid) ;
113
% determine the breakeven point where using inv(A) is faster
114
nrhs = max (1, ceil ((t3 - t1) / (t2 - t4))) ;
115
if (t1 < t3 & t2 < t4) %#ok
116
fprintf ('inv(A) breakeven: never\n') ;
117
elseif (t3 < t1 & t4 < t2) %#ok
118
fprintf ('inv(A) breakeven: > 0\n') ;
119
elseif (t1 < t3 & t4 < t2) %#ok
120
fprintf ('inv(A) breakeven: > %d\n', nrhs) ;
122
fprintf ('inv(A) breakeven: < %d\n', nrhs) ;
127
% inv(A) probably ran out of memory
128
fprintf ('inv(A) failed\n') ;
129
fprintf ('inv(A) breakeven: never\n') ;
132
%-------------------------------------------------------------------------------
134
%-------------------------------------------------------------------------------
136
% do this several times for accurate timings
141
x = A\b ; % solve Ax=b
144
trials = trials + 1 ;
148
resid = norm (A*x-b,1) / (norm (A,1) * norm (x,1) + norm (b,1)) ;
150
fprintf ('%-16s total time: %10.6f resid: %8.2e\n', ...
151
'x = A\b', t0, resid) ;
added
removed
LINFACTOR/lintest.m
