1function [err x1 x2 xset] = ltest2 (LD, L, D, L2, P, p, b, err)
2%LTEST2 test the lsubsolve mexFunction
3% Example:
4%   [err x1 x2 xset] = ltest2 (LD, L, D, L2, P, p, b, err)
5%
6% See also cholmod_test, ltest, ltest2
7
8% Copyright 2013, Timothy A. Davis, http://www.suitesparse.com
9
10if (~isreal (L))
11    b = b * (pi + 1i) ;
12end
13
14for sys = 0:8
15
16    %---------------------------------------------------------------------------
17    % solve with LDL' = A
18    %---------------------------------------------------------------------------
19
20    % solve for all of x
21    switch sys
22
23        case 0
24            x1 = P' * (L' \ (D \ (L \ (P * b)))) ;        % solve Ax = b
25
26        case 1
27            x1 =      (L' \ (D \ (L \ (    b)))) ;        % solve LDL'x = b
28
29        case 2
30            x1 =      (     (D \ (L \ (    b)))) ;        % solve LDx = b
31
32        case 3
33            x1 =      (L' \ (D \ (    (    b)))) ;        % solve DL'x = b
34
35        case 4
36            x1 =      (     (    (L \ (    b)))) ;        % solve Lx = b
37
38        case 5
39            x1 =      (L' \ (    (    (    b)))) ;        % solve L'x = b
40
41        case 6
42            x1 =      (     (D \ (    (    b)))) ;        % solve Dx = b
43
44        case 7
45            x1 =      (     (    (    (P * b)))) ;        % x = Pb
46
47        case 8
48            x1 = P' * (     (    (    (    b)))) ;        % x = P'b
49    end
50
51    % solve only for entries in xset, using lsubsolve.
52    % xset contains the pattern of b, and the reach of b in the graph of L
53    kind = 1 ;  % LDL'
54    [x2 xset] = lsubsolve (LD, kind, p, b, sys) ;
55    xset = xset'' ;
56    spok (xset) ;
57    err = max (err, norm ((x1 - x2).*xset, 1) / norm (x1,1)) ;
58    if (err > 1e-12)
59        sys
60        warning ('LDL''!') ;
61        return
62    end
63
64    %---------------------------------------------------------------------------
65    % solve with L2*L2' = A
66    %---------------------------------------------------------------------------
67
68    % solve for all of x
69    switch sys
70
71        case 0
72            x1 = P' * (L2' \ (    (L2 \ (P * b)))) ;        % solve Ax = b
73
74        case 1
75            x1 =      (L2' \ (    (L2 \ (    b)))) ;        % solve L2L2'x = b
76
77        case 2
78            x1 =      (      (    (L2 \ (    b)))) ;        % solve L2x = b
79
80        case 3
81            x1 =      (L2' \ (    (     (    b)))) ;        % solve L2'x = b
82
83        case 4
84            x1 =      (      (    (L2 \ (    b)))) ;        % solve L2x = b
85
86        case 5
87            x1 =      (L2' \ (    (     (    b)))) ;        % solve L2'x = b
88
89        case 6
90            x1 =      (      (    (     (    b)))) ;        % solve Dx = b
91
92        case 7
93            x1 =      (      (    (     (P * b)))) ;        % x = Pb
94
95        case 8
96            x1 = P' * (      (    (     (    b)))) ;        % x = P'b
97    end
98
99    % solve only for entries in xset, using lsubsolve.
100    % xset contains the pattern of b, and the reach of b in the graph of L2
101    kind = 0 ;  % L2*L2'
102    [x2 xset] = lsubsolve (L2, kind, p, b, sys) ;
103    xset = xset'' ;
104    spok (xset) ;
105    err = max (err, norm ((x1 - x2).*xset, 1) / norm (x1,1)) ;
106    if (err > 1e-12)
107        sys
108        warning ('LL''!') ;
109        return
110    end
111
112end
113