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