1## Copyright 1990-2000 Institut für Angewandte Mathematik,
2##                     Universität Karlsruhe, Germany
3## Copyright 2000-2014 Wissenschaftliches Rechnen/Softwaretechnologie,
4##                     Universität Wuppertal, Germany
5## Copyright 2015-2016 Oliver Heimlich
6##
7## This program is derived from FastLSS in CXSC, C++ library for eXtended
8## Scientific Computing (V 2.5.4), which is distributed under the terms of
9## LGPLv2+.  Original Author is Michael Zimmer.  Migration to Octave code has
10## been performed by Oliver Heimlich.
11##
12## This program is free software; you can redistribute it and/or modify
13## it under the terms of the GNU General Public License as published by
14## the Free Software Foundation; either version 3 of the License, or
15## (at your option) any later version.
16##
17## This program is distributed in the hope that it will be useful,
18## but WITHOUT ANY WARRANTY; without even the implied warranty of
19## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
20## GNU General Public License for more details.
21##
22## You should have received a copy of the GNU General Public License
23## along with this program; if not, see <http://www.gnu.org/licenses/>.
24
25## -*- texinfo -*-
26## @documentencoding UTF-8
27## @defop Method {@@infsup} mldivide (@var{X}, @var{Y})
28## @defopx Operator {@@infsup} {@var{X} \ @var{Y}}
29##
30## Return the interval matrix left division of @var{X} and @var{Y}.
31##
32## Accuracy: The result is a valid enclosure.
33##
34## @example
35## @group
36## infsup ([1, 0; 0, 2]) \ [2, 0; 0, 4]
37##   @result{} ans = 2×2 interval matrix
38##      [2]   [0]
39##      [0]   [2]
40## @end group
41## @end example
42## @seealso{@@infsup/mtimes, @@infsup/gauss}
43## @end defop
44
45## Author: Oliver Heimlich
46## Keywords: interval
47## Created: 2015-02-17
48
49function result = mldivide (A, b)
50
51  if (nargin ~= 2)
52    print_usage ();
53    return
54  endif
55  if (not (isa (A, "infsup")))
56    A = infsup (A);
57  endif
58  if (not (isa (b, "infsup")))
59    b = infsup (b);
60  elseif (isa (b, "infsupdec"))
61    ## Workaround for bug #42735
62    result = mldivide (A, b);
63    return
64  endif
65
66  if (isscalar (A.inf))
67    result = ldivide (A, b);
68    return
69  elseif (not (issquare (A.inf)))
70    error ("interval:InvalidOperand", "mldivide: Matrix is not square");
71  elseif (rows (A.inf) ~= rows (b.inf))
72    error ("interval:InvalidOperand", ["mldivide: ", ...
73                                       "nonconformant arguments ", ...
74                                       "(op1 is " num2str(rows (A.inf)) "×", ...
75                                       num2str(columns (A.inf)) ",", ...
76                                       " op2 is " num2str(rows (b.inf)) "×", ...
77                                       num2str(columns (b.inf)) ")"]);
78  elseif (isempty (A.inf))
79    result = infsup (zeros (0, columns (b.inf)));
80    return
81  endif
82
83  ## Maximum number of iterations during the verification step
84  cfg.maxIterVer = 5;
85  ## Epsilon for the verification step
86  ## (used for the Epsilon inflation during verification)
87  cfg.epsVer = 1e-1;
88  ## Maximum number of iterations during the refinement step
89  cfg.maxIterRef = 5;
90  ## Epsilon for the refinement step (stopping criterion for refinement step)
91  cfg.epsRef = 1e-5;
92  ## Maximum number of iterations during residual correction
93  ## (not available for K=1)
94  cfg.maxIterResCorr = 10;
95  ## Epsilon for the residual correction
96  cfg.epsResCorr = 1e-6;
97
98  ## An approximate inverse R of A is computed.  Then an approximate
99  ## solution x1 is computed applying a conventional residual
100  ## iteration.  For the final verification, an interval residual
101  ## iteration is performed.  An enclosure of the unique solution is
102  ## returned.
103  ##
104  ## If this first step fails, the solver will try to compute a
105  ## verified solution by using an approximate inverse of double
106  ## length.  This second step takes considerably longer than the
107  ## first step, because all computations must be performed using high
108  ## precision scalar products.
109
110  ## Compute midpoints of A and b for future reference
111  Am = mid (A);
112  bm = mid (b);
113
114  ## Approximate inversion (non-interval computation)
115  [R, cond] = inv (Am);
116  if (cond == 0)
117    result = gauss (A, b);
118    return
119  endif
120
121  ## Part 1 ====================================================================
122
123  ## Approximate solution x1 (non-interval computation)
124  x1 = R * bm;
125  x1 += R * (bm - Am * x1);
126
127  ## Interval residual x
128  x = mtimes (R, b - mtimes (A, x1, "valid"), "valid");
129
130  C = eye (rows (R)) - mtimes (R, A, "valid");
131
132  ## Verify solution x1 + x
133  [x, verified] = verify_and_refine (x, C, cfg, "valid");
134  if (verified)
135    result = x1 + x;
136    return
137  endif
138
139  ## Part 2 ====================================================================
140
141  ## R2 = inv (R * Am), with correctly rounded dot product
142  R2 = zeros (size (R));
143  for i = 1 : rows (R)
144    for j = 1 : columns (R)
145      R2 (i, j) = mpfr_vector_dot_d (0.5, R (i, :), Am (:, j));
146    endfor
147  endfor
148  [R2, cond] = inv (R2);
149  if (cond == 0)
150    result = gauss (A, b);
151    return
152  endif
153
154  ## R = R2 * R with correctly rounded dot product; error in R2
155  R1_ = R2_ = zeros (size (R));
156  for i = 1 : rows (R)
157    for j = 1 : columns (R)
158      [R1_(i, j), R2_(i, j)] = mpfr_vector_dot_d (0.5, R2 (i, :), R (:, j));
159    endfor
160  endfor
161  R = R1_; R2 = R2_; clear R1_ R2_;
162
163  ## Loop over all right hand sides
164  C_computed = false ();
165  result = infsup (zeros (size (b.inf)));
166  for s = 1 : columns (b.inf)
167    s_idx.type = "()";
168    s_idx.subs = {":", s};
169
170    ## x1 = R * bm + R2 * bm with correctly rounded dot product; error in x0
171    x1 = x0 = zeros (rows (R), 1);
172    parfor i = 1 : rows (R)
173      [x1(i), x0(i)] = mpfr_vector_dot_d (0.5, [R(i, :), R2(i, :)], ...
174                                          [bm(:, s); bm(:, s)]);
175    endparfor
176
177    ## Residual iteration (non-interval computation)
178    for k = 1 : cfg.maxIterResCorr
179      ## d = bm - Am * x1 - Am * x0 with correctly rounded dot product
180      d = zeros (rows (R), 1);
181      parfor i = 1 : rows (R)
182        d (i) = mpfr_vector_dot_d (0.5, ...
183                                   [bm(i, s), Am(i, :), Am(i, :)], ...
184                                   [1;        -x1;      -x0]);
185      endparfor
186
187      ## y0 = x0 + R * d + R2 * d with correctly rounded dot product
188      y0 = zeros (rows (R), 1);
189      parfor i = 1 : rows (R)
190        y0 (i) = mpfr_vector_dot_d (0.5, ...
191                                    [x0(i), R(i, :), R2(i, :)], ...
192                                    [1;     d;       d]);
193      endparfor
194
195      d = x1 + y0;
196      p = relative_error (d, x1 + x0);
197
198      if (p >= cfg.epsResCorr && k < cfg.maxIterResCorr)
199        ## x0 = x1 + x0 - d with correctly rounded sum
200        parfor i = 1 : rows (R)
201          x0 (i) = mpfr_vector_sum_d (0.5, [x1(i), x0(i), -d(i)]);
202        endparfor
203      endif
204
205      x1 = d;
206
207      if (p < cfg.epsResCorr)
208        break
209      endif
210    endfor
211
212    ## compute enclosure y+Y1 of the residuum b-A*x1 of the approximation x1
213    ## and initialize x:= (R+R2)*(b-A*x1), C:= I-(R+R2)*A
214
215    ## y = mid (b - A * x1)
216    y = mid ([subsref(b, s_idx), A] * [1; -x1]);
217
218    ## Y1 = b - A * x1 - y
219    Y1 = [subsref(b, s_idx), A, y] * [1; -x1; -1];
220
221    ## x = R * y + R2 * y + R * Y1 + R2 * Y1
222    x = [R, R2, R, R2] * [y; y; Y1; Y1];
223
224    ## Verifying solution x1 + x ...
225    if (all (x.inf == 0 & x.sup == 0))
226      ## exact solution! (however, not necessarily unique!)
227      subsasgn (result, s_idx, x1);
228      continue
229    endif
230
231    if (not (C_computed))
232      ## C = I - R * A - R2 * A (lazy computation)
233      C = [eye(rows (R)), R, R2] * [eye(rows (R)); -A; -A];
234      C_computed = true ();
235    endif
236
237    [x, verified] = verify_and_refine (x, C, cfg, "tight");
238    if (not (verified))
239      error ("Verification failed")
240    endif
241
242    ## The exact solution lies x1 + x
243    subsasgn (result, s_idx, x1 + x);
244  endfor
245
246endfunction
247
248## Perform an epsilon inflation
249function y = blow (x, eps)
250  y = nextout ((1 + eps) .* x - eps .* x);
251endfunction
252
253## Compute component-wise the maximum relative error
254function e = relative_error (new, old)
255  nonzero = old ~= 0 & (1e6 * abs (new) >= abs (old));
256  e = max (abs ((new (nonzero) - old (nonzero)) ./ old (nonzero)));
257
258  if (isempty (e))
259    e = 0;
260  endif
261endfunction
262
263## Interval iteration until inclusion is obtained (or max. iteration count)
264function [x, verified] = verify_and_refine (x0, C, cfg, accuracy)
265  verified = false ();
266  x = x0;
267  for p = 1 : cfg.maxIterVer
268    y = blow (x, cfg.epsVer); # epsilon inflation
269    x = x0 + mtimes (C, y, accuracy); # new iterate
270
271    verified = all (all (subset (x, y)));
272    if (verified)
273      break
274    endif
275  endfor
276
277  if (verified)
278    ## Iterative refinement
279    for p = 1 : cfg.maxIterRef
280      y = x;
281      x = intersect (x0 + mtimes (C, x, accuracy), x);
282
283      if (p == cfg.maxIterRef)
284        break
285      endif
286
287      distance = max (abs (x.inf - y.inf), ...
288                      abs (x.sup - y.sup));
289      if (max (max (distance)) <= cfg.epsRef)
290        break
291      endif
292    endfor
293  endif
294endfunction
295
296%!# unique solution
297%!assert (infsup ([1, 0; 0, 2]) \ [2, 0; 0, 4] == [2, 0; 0 2]);
298%!# no solution
299%!assert (all (isempty (infsup ([1, 0; 2, 0]) \ [3; 0])));
300%!# many solutions
301%!assert (infsup ([1, 0; 2, 0]) \ [4; 8] == infsup ([4; -inf], [4; inf]));
302%!assert (all (subset (infsup ([2, -1; -1, 2], [4, 1; 1, 4]) \ infsup ([-3; .8], [3; .8]), infsup ([-2.3; -1.1], [2.3; 1.6]))));
303