1function test9 2%TEST9 test cs_qr 3% 4% Example: 5% test9 6% See also: testall 7 8% Copyright 2006-2012, Timothy A. Davis, http://www.suitesparse.com 9 10rand ('state', 0) ; 11 12index = ssget ; 13[ignore f] = sort (max (index.nrows, index.ncols)) ; 14f = f (1:100) ; 15clf 16 17% f = 185 ; 18% f = 449 ; 19% f = 186 20 21for i = f 22 Prob = ssget (i) ; 23 disp (Prob) ; 24 A = Prob.A ; 25 [m n] = size (A) ; 26 if (m < n) 27 A = A' ; 28 end 29 [m n] = size (A) ; 30 sp = sprank (A) ; 31% if (sprank (A) < min (m,n)) 32% continue 33% end 34 35 36 for cmplex = 0:double(~ispc) 37 38 if (cmplex) 39 A = A + 1i * sprand (A) ; 40 end 41 42 Aorig = A ; 43 44 A = A (:, colamd (A)) ; 45 s1 = svd (full (A)) ; 46 47 if (cmplex) 48 try 49 tic ; 50 R = chol (A'*A) ; 51 t1 = toc ; %#ok 52 catch 53 fprintf ('chol (A''*A) failed\n') ; 54 R = [ ] ; 55 end 56 else 57 tic ; 58 R = qr (A) ; 59 t1 = toc ; %#ok 60 end 61 62 % tic ; 63 % [Q,R] = qr (A) ; 64 % t1 = toc ; 65 66 [c,h,parent] = symbfact (A, 'col') ; 67 rnz = sum (c) ; %#ok 68 tic ; 69 [V2,Beta2,p,R2] = cs_qr (sparse(A)) ; 70 t2 = toc ; %#ok 71 72 v2 = full (V2) ; 73 if (any (spones (v2) ~= spones (V2))) 74 error ('got zeros!') ; 75 end 76 77 C = A ; 78 m2 = size (V2,1) ; 79 if (m2 > m) 80 C = [A ; sparse(m2-m, n)] ; 81 end 82 C = C (p,:) ; 83 84 % [H1,R1] = myqr (C) ; 85 % err1 = norm (R1-R2,1) / norm (R1) 86 % % [svd(A) svd(R1) svd(full(R2))] 87 % s2 = svd (full (R2)) ; 88 % err2 = norm (s1 - s2) / s1 (1) 89 % fprintf ('%10.6f %10.6f cs speedup %8.3f sprank %d n %d\n', ... 90 % t1, t2, t1/t2, sp, n) ; 91 % err2 92 93 % left-looking: 94 [V,Beta3,R3] = qr_left (C) ; %#ok 95 s3 = svd (full (R2)) ; 96 err3 = norm (s1 - s3) / s1 (1) ; 97 disp ('err3 = ') ; disp (err3) ; 98 if (err3 > 1e-12) 99 error ('!') ; 100 end 101 102 % right-looking: 103 [V,Beta4,R4] = qr_right (C) ; %#ok 104 s4 = svd (full (R2)) ; 105 err4 = norm (s1 - s4) / s1 (1) ; 106 disp ('err4 = ') ; disp (err4) ; 107 if (err4 > 1e-12) 108 error ('!') ; 109 end 110 111 % H2 = full (H2) 112 % R2 = full (R2) 113 114 subplot (2,4,1) ; spy (A) ; title ('A colamd') ; 115 subplot (2,4,2) ; spy (C) ; title ('A rperm') ; 116 subplot (2,4,3) ; treeplot (parent) ; 117 subplot (2,4,4) ; spy (Aorig) ; title ('Aorig') ; 118 subplot (2,4,5) ; spy (abs(R2)>0) ; title ('spqr R, no zeros') ; 119 subplot (2,4,6) ; spy (R) ; title ('matlab R') ; 120 subplot (2,4,7) ; spy (R2) ; title ('spqr R') ; 121 subplot (2,4,8) ; spy (V2) ; title ('spqr V') ; 122 drawnow 123 124 % if (err2 > 1e-9) 125 % error ('!') ; 126 % end 127 if (m2 > m) 128 fprintf ('added %d rows, sprank %d n %d\n', m2-m, sp, n) ; 129 end 130 % pause 131 end 132 133end 134