1%% SLIP_DEMO a demo of SLIP_backslash 2% SLIP_LU is a package for solving sparse linear systems of equations 3% with a roundoff-free integer-preserving method. The result is 4% always exact, unless the matrix A is perfectly singular. 5% 6% See also vpa, SLIP_backslash, SLIP_install, SLIP_test. 7% 8% SLIP_LU: (c) 2019-2020, Chris Lourenco, Jinhao Chen, Erick Moreno-Centeno, 9% Timothy A. Davis, Texas A&M University. All Rights Reserved. See 10% SLIP_LU/License for the license. 11 12format compact 13 14%% SLIP_backslash vs MATLAB backslash: first example 15% In this first example, x = SLIP_backslash (A,b) returns an approximate 16% solution, but not because it was computed incorrectly in SLIP_backslash. 17% It is computed exactly as a rational result in SLIP_backslash with 18% arbitrary precision, but then converted to double precision on output. 19 20format long g 21load west0479 22A = west0479 ; 23n = size (A, 1) ; 24xtrue = rand (n,1) ; 25b = A*xtrue ; 26x = SLIP_backslash (A, b) ; 27% error is nonzero: x is computed exactly in rational arbitrary-precision, 28% but then lost precision when returned to MATLAB: 29err_slip = norm (x-xtrue) 30x = A\b ; 31% error is nonzero: MATLAB x=A\b experiences floating-point error 32% throughout its computations: 33err_matlab = norm (x-xtrue) 34 35%% SLIP_backslash: exact, vs MATLAB backslash: approximate 36% In this example, x = SLIP_backslash (A,b) is returned exactly in the 37% MATLAB vector x, because x contains only integers representable exactly 38% in double precision. x = A\b results in floating-point roundoff error. 39 40amax = max (abs (A), [ ], 'all') ; 41A = floor (2^20 * (A / amax)) + n * speye (n) ; 42xtrue = floor (64 * xtrue) ; 43b = A*xtrue ; 44x = SLIP_backslash (A, b) ; 45% error will be exactly zero: 46err_slip = norm (x-xtrue) 47x = A\b ; 48% error will be small but nonzero: 49err_matlab = norm (x-xtrue) 50 51%% SLIP_backslash on ill-conditioned problems 52% x = SLIP_backslash (A,b) is able to solve problems that x=A\b cannot. 53% Consider the following matrix in the MATLAB gallery: 54 55[U, b] = gallery ('wilk', 3) 56 57%% vpa can find a good but not perfect solution: 58xvpa = vpa (U) \ b 59 60% but MATLAB's numerical x = U\b computes a poor solution: 61xapprox = U \ b 62 63%% SLIP_backslash computes the exact answer 64% It returns it to MATLAB as a double vector, obtaining the exact results, 65% except for a final floating-point error in x(2): 66 67xslip = SLIP_backslash (U, b) 68err = xvpa - xslip 69relerr = double (err (2:3) ./ xvpa (2:3)) 70 71%% SLIP_backslash with exact results 72% SLIP_backslash can also return x as a cell array of strings, which 73% preserves the exact rational result. The printing option is also 74% enabled in this example. The floating-point matrices U and b are 75% converted into a scaled integer matrix before solving U*x=b with 76% SLIP LU. 77% 78% The value U(1,2)=0.9 is a floating-point number, and 0.9 cannot be 79% exactly represented in IEEE floating-point representation. It is 80% converted exactly into the rational number, 81% fl(0.9) = 45000000000000001 / 50000000000000000. 82 83option.print = 3 ; % also print the details 84option.solution = 'char' ; % return x as a cell array of strings 85 86%% 87 88xslip = SLIP_backslash (U, b, option) 89 90%% Converting an exact rational result to vpa or double 91% If SLIP_backslash returns x as a cell array of strings, it cannot 92% be immediately used in computations in MATLAB. It can be converted 93% into a vpa or double matrix, as illustrated below. The solution 94% differs slightly from the vpa solution xvpa = vpa (U)\b, since 95% the MATLAB vpa converts fl(0.9) into a decimal representation 0.9, 96% or exactly 9/10; this is not exactly equal to fl(0.9), since the 97% value 9/10 is not representable in IEEE floating-point. SLIP_backslash, 98% by contrast, converts fl(0.9) into its exact rational representation, 99% 45000000000000001 / 50000000000000000. 100 101xslip_as_vpa = vpa (xslip) 102xslip_as_double = double (vpa (xslip)) 103xvpa_as_double = double (xvpa) 104 105%% Comparing the VPA and SLIP_BACKSLASH solutions in double 106% Both vpa(U)\b and SLIP_backslash(U,b) compute the same result 107% in the end, when their results are converted to double. 108err = xvpa_as_double - xslip_as_double 109 110