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