1function gee_its_simple_test
2%GEE_ITS_SIMPLE_TEST tests the "Gee! It's Simple!" package
3% Exhaustive test of the "Gee! It's Simple!" package.  Returns the largest
4% relative residual for any solution to A*x=b (using the inf norm).  This test
5% exercises all statements in the package.  Note that the rand state is
6% modified.
7%
8% Example:
9%   gee_its_simple_test ;
10%
11% See also: gee_its_simple, gee_its_short, rand, mldivide, gee_its_simple_resid
12
13% Copyright 2007, Timothy A. Davis, http://www.suitesparse.com
14
15%-------------------------------------------------------------------------------
16% error-handling tests
17%-------------------------------------------------------------------------------
18
19fprintf ('\nTesting error handling (expect error and warning messages):\n\n');
20
21gunk = 0 ;
22ok = 0 ;
23
24lasterr ('') ;
25
26try
27    % too many inputs
28    gee_its_simple_factorize (A,gunk) ;                                     %#ok
29catch
30    ok = ok + 1 ;
31    disp (lasterr) ;
32end
33fprintf ('\n') ;
34
35try
36    % too many outputs
37    [LU,p,rcnd,gunk] = gee_its_simple_factorize (A) ;                       %#ok
38catch
39    ok = ok + 1 ;
40    disp (lasterr) ;
41end
42fprintf ('\n') ;
43
44try
45    % too few inputs
46    gee_its_simple_factorize ;                                              %#ok
47catch
48    ok = ok + 1 ;
49    disp (lasterr) ;
50end
51fprintf ('\n') ;
52
53try
54    % too many inputs
55    x = gee_its_simple (A,b,gunk) ;                                         %#ok
56catch
57    ok = ok + 1 ;
58    disp (lasterr) ;
59end
60fprintf ('\n') ;
61
62try
63    % too many outputs
64    [x,rcnd,gunk] = gee_its_simple (A,b) ;                                  %#ok
65catch
66    ok = ok + 1 ;
67    disp (lasterr) ;
68end
69fprintf ('\n') ;
70
71try
72    % too few inputs
73    x = gee_its_simple ;                                                    %#ok
74catch
75    ok = ok + 1 ;
76    disp (lasterr) ;
77end
78fprintf ('\n') ;
79
80try
81    % too many inputs
82    x = gee_its_simple_forwardsolve (A,b,gunk) ;                            %#ok
83catch
84    ok = ok + 1 ;
85    disp (lasterr) ;
86end
87fprintf ('\n') ;
88
89try
90    % too many outputs
91    [x,gunk] = gee_its_simple_forwardsolve (A,b) ;                          %#ok
92catch
93    ok = ok + 1 ;
94    disp (lasterr) ;
95end
96fprintf ('\n') ;
97
98try
99    % too few inputs
100    x = gee_its_simple_forwardsolve ;                                       %#ok
101catch
102    ok = ok + 1 ;
103    disp (lasterr) ;
104end
105fprintf ('\n') ;
106
107try
108    % too many inputs
109    x = gee_its_simple_backsolve (A,b,gunk) ;                               %#ok
110catch
111    ok = ok + 1 ;
112    disp (lasterr) ;
113end
114fprintf ('\n') ;
115
116try
117    % too many outputs
118    [x,gunk] = gee_its_simple_backsolve (A,b) ;                             %#ok
119catch
120    ok = ok + 1 ;
121    disp (lasterr) ;
122end
123fprintf ('\n') ;
124
125try
126    % too few inputs
127    x = gee_its_simple_backsolve ;                                          %#ok
128catch
129    ok = ok + 1 ;
130    disp (lasterr) ;
131end
132fprintf ('\n') ;
133
134try
135    % rectangular A
136    x = gee_its_simple (eye (4,3), ones (4,1)) ;                            %#ok
137catch
138    ok = ok + 1 ;
139    disp (lasterr) ;
140end
141fprintf ('\n') ;
142
143try
144    % A is 3D
145    x = gee_its_simple (ones (9,3,3), ones (9,1)) ;                         %#ok
146catch
147    ok = ok + 1 ;
148    disp (lasterr) ;
149end
150fprintf ('\n') ;
151
152try
153    % b is 3D
154    x = gee_its_simple (eye (3,3), ones (3,3,3)) ;                          %#ok
155catch
156    ok = ok + 1 ;
157    disp (lasterr) ;
158end
159fprintf ('\n') ;
160
161try
162    % dimensions of A and b do not matrix
163    x = gee_its_simple (eye (3,3), ones (4,1)) ;                            %#ok
164catch
165    ok = ok + 1 ;
166    disp (lasterr) ;
167end
168fprintf ('\n') ;
169
170try
171    % dimensions of L and b do not matrix
172    x = gee_its_simple_forwardsolve (eye (3,3), ones (4,1)) ;               %#ok
173catch
174    ok = ok + 1 ;
175    disp (lasterr) ;
176end
177fprintf ('\n') ;
178
179try
180    % dimensions of U and b do not matrix
181    x = gee_its_simple_backsolve (eye (3,3), ones (4,1)) ;                  %#ok
182catch
183    ok = ok + 1 ;
184    disp (lasterr) ;
185end
186fprintf ('\n') ;
187
188% singular matrix
189lastwarn ('') ;
190x = gee_its_simple (0, 1) ;                                                 %#ok
191[msg, id] = lastwarn ;
192if (~isempty (msg) & ~isempty (id))                                         %#ok
193    ok = ok + 1 ;
194end
195fprintf ('\n') ;
196
197% ill-conditioned matrix
198lastwarn ('') ;
199x = gee_its_simple ([1e30 2e30 ; 1 1], [1 ; 1]) ;                           %#ok
200[msg, id] = lastwarn ;
201if (~isempty (msg) & ~isempty (id))                                         %#ok
202    ok = ok + 1 ;
203end
204
205if (ok ~= 20)
206    error ('test failed') ;
207end
208
209fprintf ('\n\nError-handing tests complete (all error messages and warnings\n');
210fprintf ('shown above were expected).  Now testing for accuracy:\n\n') ;
211
212%-------------------------------------------------------------------------------
213% compare accuracy vs. backslash
214%-------------------------------------------------------------------------------
215
216maxerr1 = 0 ;   % largest residual for A\b (gee_its_sweet)
217maxerr2 = 0 ;   % largest residual for gee_its_simple (A,b)
218maxerr3 = 0 ;   % largest residual for gee_its_short (A,b)
219maxerr4 = 0 ;   % largest residual for gee_its_too_short (A,b)
220rmax = 0 ;      % largest relative difference in rcond
221
222nmax = 50 ;     % largest dimension of A to test
223cmax = 10 ;     % largest number of columns of b to test
224ntrials = 2 ;   % number of trials for each x=A\b
225
226rand ('state', 0) ;
227
228for n = 0:nmax
229    fprintf ('.') ;
230    for c = 0:cmax
231        for trial = 1:ntrials
232
233            % set up the system
234            A = rand (n) ;
235            b = rand (n,c) ;
236
237            % solve it four different ways
238            x1 = gee_its_sweet (A,b) ;      % this is just a one-liner: x=A\b
239            x2 = gee_its_simple (A,b) ;
240            x3 = gee_its_short (A,b) ;
241            x4 = gee_its_too_short (A,b) ;
242
243            % get the relative residuals
244            err1 = gee_its_simple_resid (A, x1, b) ;
245            err2 = gee_its_simple_resid (A, x2, b) ;
246            err3 = gee_its_simple_resid (A, x3, b) ;
247            err4 = gee_its_simple_resid (A, x4, b) ;
248
249            maxerr1 = max (maxerr1, err1) ;
250            maxerr2 = max (maxerr2, err2) ;
251            maxerr3 = max (maxerr3, err3) ;
252            maxerr4 = max (maxerr4, err4) ;
253
254            if (max ([err1 err2 err3]) > 1e-14)
255                error ('test failed') ;
256            end
257
258            % test rcond
259            if (n > 0)
260                [L,U,p] = lu (A) ;                                          %#ok
261                r1 = min (abs (diag (U))) / max (abs (diag (U))) ;
262                [LU,p,r2] = gee_its_simple_factorize (A) ;
263                if (r1 ~= 0)
264                    r = abs (r1 - r2) / r1 ;
265                    rmax = max (rmax, r) ;
266                end
267                if (r > 1e-10)
268                    error ('test failed') ;
269                end
270            end
271        end
272    end
273end
274fprintf ('\n') ;
275
276fprintf ('max residual for backslash:         %g\n', maxerr1) ;
277fprintf ('max residual for gee_its_simple:    %g\n', maxerr2) ;
278fprintf ('max residual for gee_its_short:     %g\n', maxerr3) ;
279fprintf ('max residual for gee_its_too_short: %g (no pivoting!)\n', maxerr4) ;
280
281fprintf ('\n\nAll tests passed OK\n') ;
282
283