1function test8 2%TEST8 test cs_cholsol, cs_lusol 3% 4% Example: 5% test8 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 14% f = f(1) 15 16for i = f 17 Prob = ssget (i) ; 18 disp (Prob) ; 19 A = Prob.A ; 20 [m n] = size (A) ; 21 if (m ~= n) 22 continue 23 end 24 25 for cmplex = 0:double(~ispc) 26 27 if (cmplex) 28 A = A + 0.1i * (sprand (tril (A,-1) + triu (A,1))) ; 29 end 30 31 spd = 0 ; 32 if (m == n) 33 if (nnz (A-A') == 0) 34 try 35 p = amd (A) ; 36 catch 37 p = symamd (A) ; 38 end 39 [R,p] = chol (A (p,p)) ; 40 spd = (p == 0) ; 41 end 42 end 43 44 if (spd) 45 C = A ; 46 else 47 C = A*A' + n*speye (n) ; 48 try 49 p = amd (C) ; 50 catch 51 p = symamd (C) ; 52 end 53 try 54 R = chol (C (p,p)) ; %#ok 55 catch 56 continue 57 end 58 end 59 60 b = rand (n,1) ; 61 62 x1 = C\b ; 63 x2 = cs_cholsol (C,b) ; 64 r1 = norm (C*x1-b,1) / norm (C,1) ; 65 r2 = norm (C*x2-b,1) / norm (C,1) ; 66 err = abs (r1-r2) ; 67 fprintf ('err %g\n', err) ; 68 if (err > 1e-10) 69 error ('!') ; 70 end 71 72 x2 = cs_lusol (C,b, 1, 0.001) ; 73 r2 = norm (C*x2-b,1) / norm (C,1) ; 74 err = abs (r1-r2) ; 75 fprintf ('err %g (lu with amd(A+A'')\n', err) ; 76 if (err > 1e-10) 77 error ('!') ; 78 end 79 80 if (m ~= n) 81 continue ; 82 end 83 84 x1 = A\b ; 85 r1 = norm (A*x1-b,1) / norm (A,1) ; 86 if (r1 < 1e-6) 87 x2 = cs_lusol (A,b) ; 88 r2 = norm (A*x2-b,1) / norm (A,1) ; 89 fprintf ('lu resid %g %g\n', r1, r2) ; 90 end 91 end 92end 93