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