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