1function test6
2%TEST6 test cs_reach, cs_reachr, cs_lsolve, cs_usolve
3%
4% Example:
5%   test6
6% See also: testall
7
8% Copyright 2006-2012, Timothy A. Davis, http://www.suitesparse.com
9
10rand ('state', 0)
11maxerr = 0 ;
12clf
13for trial = 1:201
14    n = fix (100 * rand (1)) ;
15    d = 0.1 * rand (1) ;
16    L = tril (sprandn (n,n,d),-1) + sprand (speye (n)) ;
17    b = sprandn (n,1,d) ;
18
19    if (~ispc)
20        if (rand ( ) > .5)
21            L = L + 1i*sprand(L) ;
22        end
23        if (rand ( ) > .5)
24            b = b + 1i*sprand(b) ;
25        end
26    end
27
28    for uplo = 0:1
29
30        if (uplo == 1)
31            % solve Ux=b instead ;
32            L = L' ;
33        end
34
35        x = L\b ;
36        sr = 1 + cs_reachr (L,b) ;
37        sz = 1 + cs_reachr (L,b) ;
38
39        check_if_same (sr,sz) ;
40
41        s2 = 1 + cs_reach (L,b) ;
42
43        try
44            if (uplo == 0)
45                x3 = cs_lsolve (L,b) ;
46            else
47                x3 = cs_usolve (L,b) ;
48            end
49        catch
50            if (isreal (L) & isreal (b))                                    %#ok
51                lasterr
52                error ('!') ;
53            end
54            % punt: sparse(L)\sparse(b) not handled by cs_lsolve or cs_usolve
55            x3 = L\b ;
56        end
57
58        % x3 is NOT returned in sorted order, so it is not a valid MATLAB
59        % sparse matrix.  It might also have explicit zeros.  Double-transpose
60        % it to sort it, and multiply one to remove explicit zeros.
61        x3 = 1 * (x3')' ;
62
63        spy ([L b x x3])
64        drawnow
65
66        s = sort (sr) ;
67        [i j xx] = find (x) ;                                               %#ok
68        [i3 j3 xx3] = find (x3) ;                                           %#ok
69
70        if (isempty (i))
71            if (~isempty (s))
72                i       %#ok
73                s       %#ok
74                error ('!') ;
75            end
76        elseif (any (s ~= i))
77            i       %#ok
78            s       %#ok
79            error ('!') ;
80        end
81
82        if (isempty (i3))
83            if (~isempty (s))
84                i3      %#ok
85                s       %#ok
86                error ('!') ;
87            end
88        elseif (any (s ~= sort (i3)))
89            s       %#ok
90            i3      %#ok
91            error ('!') ;
92        end
93
94        if (any (s2 ~= sr))
95            s2      %#ok
96            sr      %#ok
97            error ('!') ;
98        end
99
100        err = norm (x-x3,1) ;
101        if (err > 1e-12)
102            x       %#ok
103            x3      %#ok
104            uplo    %#ok
105            err     %#ok
106            error ('!')
107        end
108
109        maxerr = max (maxerr, err) ;
110
111    end
112
113    drawnow
114end
115fprintf ('maxerr = %g\n', maxerr) ;
116