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 (~isreal (A) | m ~= n)                                               %#ok
22        continue
23    end
24
25    spd = 0 ;
26    if (m == n)
27        if (nnz (A-A') == 0)
28            try
29                p = amd (A) ;
30            catch
31                p = symamd (A) ;
32            end
33            [R,p] = chol (A (p,p)) ;
34            spd = (p == 0) ;
35        end
36    end
37
38    if (spd)
39        C = A ;
40    else
41        C = A*A' + n*speye (n) ;
42        try
43            p = amd (C) ;
44        catch
45            p = symamd (C) ;
46        end
47        try
48            R = chol (C (p,p)) ;
49        catch
50            continue
51        end
52    end
53
54    b = rand (n,1) ;
55
56    x1 = C\b ;
57    x2 = cs_cholsol (C,b) ;
58    r1 = norm (C*x1-b,1) / norm (C,1) ;
59    r2 = norm (C*x2-b,1) / norm (C,1) ;
60    err = abs (r1-r2) ;
61    fprintf ('err %g\n', err) ;
62    if (err > 1e-10)
63        error ('!') ;
64    end
65
66    x2 = cs_lusol (C,b, 1, 0.001) ;
67    r2 = norm (C*x2-b,1) / norm (C,1) ;
68    err = abs (r1-r2) ;
69    fprintf ('err %g (lu with amd(A+A'')\n', err) ;
70    if (err > 1e-10)
71        error ('!') ;
72    end
73
74    if (m ~= n)
75        continue ;
76    end
77
78    x1 = A\b ;
79    r1 = norm (A*x1-b,1) / norm (A,1) ;
80    if (r1 < 1e-6)
81        x2 = cs_lusol (A,b) ;
82        r2 = norm (A*x2-b,1) / norm (A,1) ;
83        fprintf ('lu resid %g %g\n', r1, r2) ;
84    end
85end
86