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