1function test7
2%TEST7 test cs_lu
3%
4% Example:
5%   test7
6% See also: testall
7
8% Copyright 2006-2012, Timothy A. Davis, http://www.suitesparse.com
9
10index = ssget ;
11[ignore f] = sort (max (index.nrows, index.ncols)) ;
12f = f (1:100) ;
13
14clf
15
16maxerr1 = 0 ;
17maxerr2 = 0 ;
18
19for i = f
20    Prob = ssget (i) ;
21    disp (Prob) ;
22    A = Prob.A ;
23
24    [m n] = size (A) ;
25    if (m ~= n)
26        continue
27    end
28
29    for cmplex = 0:1
30
31        if (cmplex)
32            A = A + norm(A,1) * sprand (A) / 3 ;
33        end
34
35        [L,U,P] = lu (A) ;
36
37        udiag = full (diag (U)) ;
38        umin = min (abs (udiag)) ;
39        fprintf ('umin %g\n', umin) ;
40
41        if (umin > 1e-14)
42
43            [L2,U2,p] = cs_lu (A) ;
44
45            subplot (3,4,1) ; spy (A) ;
46            subplot (3,4,2) ; spy (A(p,:)) ;
47            subplot (3,4,3) ; spy (L2) ;
48            subplot (3,4,4) ; spy (U2) ;
49
50            err1 = norm (L*U-P*A,1) / norm (A,1) ;
51            err2 = norm (L2*U2-A(p,:),1) / norm (A,1) ;
52            fprintf ('err %g %g\n', err1, err2) ;
53
54            if (err1 > 1e-10 | err2 > 1e-10)                                %#ok
55                error ('!') ;
56            end
57            maxerr1 = max (maxerr1, err1) ;
58            maxerr2 = max (maxerr2, err2) ;
59
60        end
61
62        q = colamd (A) ;
63
64        [L,U,P] = lu (A (:,q)) ;
65
66        udiag = full (diag (U)) ;
67        umin = min (abs (udiag)) ;
68        fprintf ('umin %g with q\n', umin) ;
69
70        if (umin > 1e-14)
71
72            [L2,U2,p,q2] = cs_lu (A) ;
73
74            subplot (3,4,5) ; spy (A) ;
75            subplot (3,4,6) ; spy (A(p,q2)) ;
76            subplot (3,4,7) ; spy (L2) ;
77            subplot (3,4,8) ; spy (U2) ;
78
79            err1 = norm (L*U-P*A(:,q),1) / norm (A,1) ;
80            err2 = norm (L2*U2-A(p,q2),1) / norm (A,1) ;
81            fprintf ('err %g %g\n', err1, err2) ;
82
83            if (err1 > 1e-10 | err2 > 1e-10)                                %#ok
84                error ('!') ;
85            end
86            maxerr1 = max (maxerr1, err1) ;
87            maxerr2 = max (maxerr2, err2) ;
88        end
89
90
91        try
92            q = amd (A) ;
93        catch
94            q = symamd (A) ;
95        end
96
97        tol = 0.01 ;
98
99        [L,U,P] = lu (A (q,q), tol) ;
100
101        udiag = full (diag (U)) ;
102        umin = min (abs (udiag)) ;
103        fprintf ('umin %g with amd q\n', umin) ;
104
105        if (umin > 1e-14)
106
107            [L2,U2,p,q2] = cs_lu (A,tol) ;
108
109            subplot (3,4,9) ; spy (A) ;
110            subplot (3,4,10) ; spy (A(p,q2)) ;
111            subplot (3,4,11) ; spy (L2) ;
112            subplot (3,4,12) ; spy (U2) ;
113
114            err1 = norm (L*U-P*A(q,q),1) / norm (A,1) ;
115            err2 = norm (L2*U2-A(p,q2),1) / norm (A,1) ;
116            lbig = full (max (max (abs (L2)))) ;
117            fprintf ('err %g %g lbig %g\n', err1, err2, lbig) ;
118            if (lbig > 1/tol)
119                error ('L!') ;
120            end
121
122            if (err1 > 1e-10 | err2 > 1e-10)                                %#ok
123                error ('!') ;
124            end
125            maxerr1 = max (maxerr1, err1) ;
126            maxerr2 = max (maxerr2, err2) ;
127        end
128
129        drawnow
130        % pause
131
132    end
133end
134
135fprintf ('maxerr %g %g\n', maxerr1, maxerr2) ;
136