1function test10 2%TEST10 test cs_qr 3% 4% Example: 5% test10 6% See also: testall 7 8% Copyright 2006-2012, Timothy A. Davis, http://www.suitesparse.com 9 10rand ('state', 0) ; 11 12% f = 185 ; 13% f = 449 ; 14clf 15 16for trials = 1:100 17 18 m = fix (100 * rand (1)) ; 19 n = fix (100 * rand (1)) ; 20 d = 0.1 * rand (1) ; 21 A = sprandn (m, n, d) ; 22 [m n] = size (A) ; 23 if (m < n) 24 A = A' ; 25 end 26 [m n] = size (A) ; 27 sp = sprank (A) ; 28 % if (sp < n) 29 % continue ; 30 % end 31 32 for cmplex = 0:double(~ispc) 33 34 if (cmplex) 35 A = A + 1i * sprand (A) * norm (A,1) / 10 ; 36 end 37 38 Aorig = A ; 39 40 % A = A (:, colamd (A)) ; 41 42 if (cmplex) 43 tic ; 44 R = chol (A'*A + speye (n)) ; 45 t1 = toc ; 46 else 47 tic ; 48 R = qr (A) ; 49 t1 = toc ; 50 end 51 52 % tic ; 53 % [Q,R] = qr (A) ; 54 % t1 = toc ; 55 56 [c,h,parent] = symbfact (A, 'col') ; %#ok 57 rnz = sum (c) ; %#ok 58 tic ; 59 [V2,Beta2,p,R2] = cs_qr (sparse(A)) ; 60 t2 = toc ; 61 62 C = A ; 63 m2 = size (V2,1) ; 64 if (m2 > m) 65 C = [A ; sparse(m2-m, n)] ; 66 end 67 C = C (p,:) ; 68 69 [H1,R1] = myqr (C) ; 70 err1 = norm (R1-R2,1) / norm (R1) ; 71 disp ('err1 = ') ; 72 disp (err1) ; 73 % [svd(A) svd(R1) svd(full(R2))] 74 s1 = svd (full (A)) ; 75 s2 = svd (full (R2)) ; 76 if (n > 0) 77 err2 = norm (s1 - s2) / s1 (1) ; 78 disp ('err2 = ') ; 79 disp (err2) ; 80 else 81 err2 = 0 ; 82 end 83 fprintf ('%10.6f %10.6f cs speedup %8.3f sprank %d vs %d\n', t1, t2, t1/t2, sp, n) ; 84 85 % H2 = full (H2) 86 % R2 = full (R2) 87 88 subplot (2,4,1) ; spy (A) ; title ('A colamd') ; 89 subplot (2,4,4) ; spy (Aorig) ; title ('Aorig') ; 90 subplot (2,4,2) ; spy (C) ; title ('A rperm') ; 91 subplot (2,4,5) ; spy (abs(R2)>0) ; title ('spqr R, no zeros') ; 92 subplot (2,4,6) ; spy (R) ; title ('matlab R') ; 93 subplot (2,4,7) ; spy (R2) ; title ('spqr R') ; 94 subplot (2,4,8) ; spy (V2) ; title ('spqr H') ; 95 drawnow 96 97 if (err2 > 1e-9) 98 error ('!') ; 99 end 100 101 if (m2 > m) 102 fprintf ('added %d rows, sprank %d n %d\n', m2-m, sp, n) ; 103 end 104 end 105end 106