inv function failing?
Enzo Alberto Dari <[email protected]> Sun, 4 Dec 2022 17:01:29 -0300
| Newsgroups | gmane.comp.gnu.octave.general |
|---|---|
| Message-ID | <CALOFp8dso9nYQ1o0DorUEcO=uE5Nh4nSX0dCyq75UDB7nnb1fA@mail.gmail.com> |
While upgrading my debian OS from 10.x to 11.x (octave 4.4.5 to 6.2.0), one of my scripts started failing. I managed to create a simple test (attached if the mailing list doesn't block it) that reproduces the problem. The matrix is blocked, so there is an alternate method to compute its inverse by blocks that always works. This test verifies if the computed inverse times the original matrix has condition number 1 (as it should be for the identity matrix). To make things more complicated, I ran the test in 5 machines with debian 11.5 (octave 6.2.0) and it failed in 3 of them, working ok in the other 2. Perhaps something to do with new processor features? (apparently is failing in machines with "newer" processors): Intel(R) Core(TM) i5 CPU 750 @ 2.67GHz: Ok. Intel(R) Core(TM) i7-4930K CPU @ 3.40GHz: Ok. Intel(R) Core(TM) i5-1035G1 CPU @ 1.00GHz: Fails. Intel(R) Core(TM) i9-7900X CPU @ 3.30GHz: Fails. Intel(R) Core(TM) i9-9900X CPU @ 3.50GHz: Fails. Any clue about what may be happening? Any tips for further analysis? Thanks in advance, -- Enzo A. Dari Profesor Titular Instituto Balseiro <http://www.ib.edu.ar/> ---------- We are transitioning to a web based forum for community help discussions at https://octave.discourse.group/c/help
chkBlockedInverse.m
(text/x-objcsrc, 3.6 KB)
% non-singular matrix b=[7110.327, -2592.219, 631.419, -288.541, 169.250, -113.431, 82.646, -63.812, 51.448, -42.914, 36.804, -32.315, 28.958, -26.424, 24.511, -23.083, 22.050, -21.350, 20.944; -1218.551, 1508.124, -720.486, 169.250, -74.433, 42.572, -28.131, 20.364, -15.701, 12.683, -10.624, 9.164, -8.102, 7.317, -6.734, 6.304, -5.996, 5.789, -5.669; 169.250, -482.641, 674.499, -350.244, 82.646, -36.010, 20.364, -13.333, 9.592, -7.371, 5.950, -4.991, 4.320, -3.838, 3.488, -3.235, 3.056, -2.936, 2.868; -49.544, 82.646, -268.958, 399.001, -215.550, 51.448, -22.463, 12.683, -8.285, 5.950, -4.571, 3.694, -3.108, 2.702, -2.415, 2.211, -2.069, 1.976, -1.923; 20.364, -27.810, 51.448, -178.696, 275.664, -152.325, 36.804, -16.173, 9.164, -6.000, 4.320, -3.329, 2.702, -2.285, 2.000, -1.803, 1.668, -1.580, 1.531; -10.260, 12.683, -18.791, 36.804, -132.946, 211.205, -118.438, 28.958, -12.831, 7.317, -4.818, 3.488, -2.705, 2.211, -1.886, 1.668, -1.522, 1.428, -1.376; 5.950, -6.944, 9.164, -14.257, 28.958, -107.451, 174.780, -99.104, 24.511, -10.963, 6.304, -4.185, 3.056, -2.392, 1.976, -1.707, 1.531, -1.420, 1.359; -3.839, 4.320, -5.318, 7.317, -11.762, 24.511, -92.803, 153.928, -88.144, 22.050, -9.966, 5.789, -3.883, 2.868, -2.273, 1.906, -1.674, 1.531, -1.453; 2.702, -2.966, 3.488, -4.451, 6.304, -10.377, 22.050, -84.854, 143.052, -82.752, 20.944, -9.576, 5.630, -3.826, 2.868, -2.312, 1.976, -1.775, 1.668; -2.050, 2.211, -2.519, 3.056, -4.000, 5.789, -9.704, 20.944, -81.727, 139.664, -81.727, 20.944, -9.704, 5.789, -4.000, 3.056, -2.519, 2.211, -2.050; 1.668, -1.775, 1.976, -2.312, 2.868, -3.826, 5.630, -9.576, 20.944, -82.752, 143.052, -84.854, 22.050, -10.377, 6.304, -4.451, 3.488, -2.966, 2.702; -1.453, 1.531, -1.674, 1.906, -2.273, 2.868, -3.883, 5.789, -9.966, 22.050, -88.144, 153.928, -92.803, 24.511, -11.762, 7.317, -5.318, 4.320, -3.839; 1.359, -1.420, 1.531, -1.707, 1.976, -2.392, 3.056, -4.185, 6.304, -10.963, 24.511, -99.104, 174.780, -107.451, 28.958, -14.257, 9.164, -6.944, 5.950; -1.376, 1.428, -1.522, 1.668, -1.886, 2.211, -2.705, 3.488, -4.818, 7.317, -12.831, 28.958, -118.438, 211.205, -132.946, 36.804, -18.791, 12.683, -10.260; 1.531, -1.580, 1.668, -1.803, 2.000, -2.285, 2.702, -3.329, 4.320, -6.000, 9.164, -16.173, 36.804, -152.325, 275.664, -178.696, 51.448, -27.810, 20.364; -1.923, 1.976, -2.069, 2.211, -2.415, 2.702, -3.108, 3.694, -4.571, 5.950, -8.285, 12.683, -22.463, 51.448, -215.550, 399.001, -268.958, 82.646, -49.544; 2.868, -2.936, 3.056, -3.235, 3.488, -3.838, 4.320, -4.991, 5.950, -7.371, 9.592, -13.333, 20.364, -36.010, 82.646, -350.244, 674.499, -482.641, 169.250; -5.669, 5.789, -5.996, 6.304, -6.734, 7.317, -8.102, 9.164, -10.624, 12.683, -15.701, 20.364, -28.131, 42.572, -74.433, 169.250, -720.486, 1508.124, -1218.551; 20.944, -21.350, 22.050, -23.083, 24.511, -26.424, 28.958, -32.315, 36.804, -42.914, 51.448, -63.812, 82.646, -113.431, 169.250, -288.541, 631.419, -2592.219, 7110.327]; sizeb=n=size(b,1) rankb=rank(b) % Builds blocked matrix B=[b zeros(n); zeros(n) eye(n)]; % Computes inverse inv1=inv(B); % Computes invese by blocks: non-trivial block: invb=inv(b); % Build inverse inv2=[invb zeros(n); zeros(n) eye(n)]; % Both inverse matrices should be equal diffinvs=norm(inv1-inv2) % All these condition numbers should be 1 cond(inv1*B) cond(B*inv1) cond(inv2*B) cond(B*inv2) % Check also with another ordering B_=[eye(n) zeros(n); zeros(n) b]; % use inv function inv1_=inv(B_); %invb=inv(b); % Build inverse by blocks inv2_=[eye(n) zeros(n); zeros(n) invb]; % Both inverse matrices should be equal diffinvs_=norm(inv1-inv2) % All these condition numbers should be 1 cond(inv1_*B_) cond(B_*inv1_) cond(inv2_*B_) cond(B_*inv2_)