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