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