1function spqr_demo
2%SPQR_DEMO short demo of SuiteSparseQR
3%
4% Example:
5%   spqr_demo
6%
7% See also SPQR, SPQR_SOLVE, SPQR_QMULT, SPQR_MAKE, SPQR_INSTALL.
8
9% Copyright 2008, Timothy A. Davis, http://www.suitesparse.com
10
11help spqr
12help spqr_solve
13help spqr_qmult
14help spqr_demo
15
16% input ('Hit enter to start the SuiteSparseQR demo: ', 's') ;
17fprintf ('\nTesting SuiteSparseQR functions ... please wait ...\n') ;
18
19load west0479 ;
20A = west0479 ;
21rand ('state', 0) ;     %#ok
22m = size (A,1) ;
23
24clf
25
26maxerr = 0 ;
27
28for acomplex = 0:1
29
30    if (acomplex)
31        A = A + 1i * sprand (A) ;
32    end
33    anorm = norm (A,1) ;
34
35    R1 = spqr (A) ;
36    err = norm (R1'*R1 - A'*A, 1) / anorm^2 ;
37    maxerr = max (maxerr, err) ;
38
39    [Q,R,E] = spqr (A) ;
40    err = norm (Q*R-A*E, 1) / anorm ;
41    maxerr = max (maxerr, err) ;
42
43    [H,R,E] = spqr (A, struct ('Q', 'Householder')) ;
44    err = norm (spqr_qmult (H,R,1) - A*E, 1) / anorm ;
45    maxerr = max (maxerr, err) ;
46
47    if (acomplex)
48        subplot (2,5,7)  ; spy (R1)  ; title ('R, no permutation (complex)') ;
49        subplot (2,5,8)  ; spy (R)   ; title ('R with colamd (complex)') ;
50        subplot (2,5,9)  ; spy (Q)   ; title ('Q with colamd (complex)') ;
51        subplot (2,5,10) ; spy (H.H) ; title ('H with colamd (complex)') ;
52    else
53        subplot (2,5,1)  ; spy (A)   ; title ('A') ;
54        subplot (2,5,2)  ; spy (R1)  ; title ('R, no permutation') ;
55        subplot (2,5,3)  ; spy (R)   ; title ('R with colamd') ;
56        subplot (2,5,4)  ; spy (Q)   ; title ('Q with colamd') ;
57        subplot (2,5,5)  ; spy (H.H) ; title ('H with colamd') ;
58    end
59    drawnow
60
61    % test spqr_solve: real/complex, sparse/full,
62    % single/multiple right-hand-sides
63    for bcomplex = 0:1
64        for bsparse = 0:1
65            for nrhs = 0:10
66                if (bsparse)
67                    b = sprand (m, nrhs, 0.1) ;
68                else
69                    b = rand (m, nrhs) ;
70                end
71                if (bcomplex)
72                    b = b + 1i * sprand (b) ;
73                end
74                x = spqr_solve (A,b) ;
75                err = norm (A*x-b,1) / max (anorm * norm(x,1) + norm (b,1), 1) ;
76                maxerr = max (maxerr, err) ;
77            end
78        end
79    end
80end
81
82fprintf ('\nQR maximum error: %g\n\n', maxerr) ;
83if (maxerr > 1e-12)
84    error ('One or more tests failed; error is high!') ;
85end
86
87% ------------------------------------------------------------------------------
88% compare spqr_solve
89% ------------------------------------------------------------------------------
90
91fprintf ('Compare performance with MATLAB on a dense least-squares problem:\n');
92A = rand (2000,1000) ;
93b = rand (2000,1) ;
94S = sparse (A) ;
95fprintf ('\nA = rand (2000,1000) ;\nb = rand (2000,1) ;\nS = sparse (A) ;\n') ;
96
97fprintf ('tic, x = spqr_solve(S,b) ; toc  ') ;
98tic
99x = spqr_solve (S, b) ;
100t1 = toc ;
101r1 = norm (A*x-b,1) ;
102r1b = norm (A'*(A*x)-A'*b,1) / norm (A'*A,1) ;
103fprintf ('%% time %8.3f residual %8.3e normal eqn resid %8.3e\n', t1,r1, r1b) ;
104
105fprintf ('tic, x = A\\b ; toc              ') ;
106tic
107x = A\b ;
108t2 = toc ;
109r2 = norm (A*x-b,1) ;
110r2b = norm (A'*(A*x)-A'*b,1) / norm (A'*A,1) ;
111fprintf ('%% time %8.3f residual %8.3e normal eqn resid %8.3e\n', t2,r2,r2b) ;
112
113fprintf ('tic, x = S\\b ; toc              ') ;
114tic
115x = S\b ;
116t3 = toc ;
117r3 = norm (A*x-b,1) ;
118r3b = norm (A'*(A*x)-A'*b,1) / norm (A'*A,1) ;
119fprintf ('%% time %8.3f residual %8.3e normal eqn resid %8.3e\n', t3,r3,r3b) ;
120
121% ------------------------------------------------------------------------------
122% compare spqr with 1000-by-2000 system
123% ------------------------------------------------------------------------------
124
125fprintf ('\nA = rand (1000,2000) ;\nS = sparse (A) ;\n') ;
126A = rand (1000,2000) ;
127S = sparse (A) ;
128
129fprintf ('tic, R = spqr(S) ; toc          ') ;
130tic
131R = spqr (S) ;
132t1 = toc ;
133r1 = norm (R'*R-A'*A,1) / norm(A,1)^2 ;
134fprintf ('%% time %8.3f error %8.3e\n', t1,r1) ;
135clear R
136
137fprintf ('tic, R = qr(A)   ; toc          ') ;
138tic
139R = qr (A) ;
140t2 = toc ;
141R = triu (R) ;
142r2 = norm (R'*R-A'*A,1) / norm(A,1)^2 ;
143fprintf ('%% time %8.3f error %8.3e\n', t2,r2) ;
144clear R
145
146fprintf ('tic, R = qr(S)   ; toc          ') ;
147tic
148R = qr (S) ;
149t3 = toc ;
150r3 = norm (R'*R-A'*A,1) / norm(A,1)^2 ;
151fprintf ('%% time %8.3f error %8.3e\n', t3,r3) ;
152clear R
153
154% ------------------------------------------------------------------------------
155% compare spqr with 100-by-20000 system
156% ------------------------------------------------------------------------------
157
158A = rand (100,20000) ;
159S = sparse (A) ;
160fprintf ('\nA = rand (100,20000) ;\nS = sparse (A) ;\n') ;
161
162fprintf ('tic, R = spqr(S) ; toc          ') ;
163tic
164R = spqr (S) ;                                                              %#ok
165t1 = toc ;
166fprintf ('%% time %8.3f\n', t1) ;
167clear R
168
169fprintf ('tic, R = qr(A)   ; toc          ') ;
170tic
171R = qr (A) ;                                                                %#ok
172t2 = toc ;
173fprintf ('%% time %8.3f\n', t2) ;
174clear R
175
176% skip the old MATLAB QR ... it's way to slow ...
177%   fprintf ('tic, R = qr(S)   ; toc          ') ;
178%   try
179%       tic
180%       R = qr (S) ;                                                        %#ok
181%       t3 = toc ;
182%       fprintf ('%% time %8.3f\n', t3) ;
183%   catch                                                                   %#ok
184%       fprintf ('%% MATLAB sparse qr failed ...\n') ;
185%   end
186%   clear R
187
188fprintf ('All spqr tests passed\n') ;
189