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