Mercurial > octave

changeset 28152:4609d001daee stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
condest.m: Fix estimate when matrix is not symmetric (bug #57968). * condest.m (solve_sparse, solve_not_sparse): In switch statement, swap the calculation method between "notransp" and "transp" cases. Add regression BIST tests.
author Rik <rik@octave.org>
date Wed, 11 Mar 2020 13:11:38 -0700
parents 026bff6a54d7
children f67f582da447 e7fe6703a81f
files scripts/linear-algebra/condest.m
diffstat 1 files changed, 22 insertions(+), 5 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/linear-algebra/condest.m	Tue Mar 10 07:51:20 2020 -0400
+++ b/scripts/linear-algebra/condest.m	Wed Mar 11 13:11:38 2020 -0700
@@ -208,7 +208,7 @@
   warning ("off", "Octave:nearly-singular-matrix", "local");
 
   if (! have_solve_normest1)
-     ## prepare solve in normest1 form
+    ## prepare solve in normest1 form
     if (issparse (A))
       [L, U, P, Pc] = lu (A);
       solve = @(flag, x) solve_sparse (flag, x, n, real_op, L, U, P, Pc);
@@ -240,15 +240,18 @@
 endfunction
 
 function value = solve_sparse (flag, x, n, real_op, L , U , P , Pc)
+  ## FIXME: Sparse algorithm is less accurate than full matrix version
+  ##        See BIST test for non-orthogonal matrix where relative tolerance
+  ##        of 1e-12 is used for sparse, but 4e-16 for full matrix.
   switch (flag)
     case "dim"
       value = n;
     case "real"
       value = real_op;
     case "notransp"
-      value = P' * (L' \ (U' \ (Pc * x)));
+      value = Pc' * (U \ (L \ (P * x)));
     case "transp"
-      value = Pc' * (U \ (L \ (P * x)));
+      value = P' * (L' \ (U' \ (Pc * x)));
   endswitch
 endfunction
 
@@ -259,9 +262,9 @@
     case "real"
       value = real_op;
     case "notransp"
-      value = P' * (L' \ (U' \ x));
+      value = U \ (L \ (P * x));
     case "transp"
-      value = U \ (L \ (P * x));
+      value = P' * (L' \ (U' \ x));
   endswitch
 endfunction
 
@@ -339,6 +342,7 @@
 %! cA_test = norm (inv (A^2), 1) * norm (A^2, 1);
 %! assert (cA, cA_test, -2^-6);
 
+## Test singular matrices
 %!test <*46737>
 %! A = [ 0         0         0
 %!       0   3.33333 0.0833333
@@ -347,6 +351,19 @@
 %! assert (cest, Inf);
 %! assert (v, []);
 
+## Test non-orthogonal matrices
+%!test <*57968>
+%! A = reshape (sqrt (0:15), 4, 4);
+%! cexp = norm (A, 1) * norm (inv (A), 1);
+%! cest = condest (A);
+%! assert (cest, cexp, -2*eps);
+
+%!test <*57968>
+%! As = sparse (reshape (sqrt (0:15), 4, 4));
+%! cexp = norm (As, 1) * norm (inv (As), 1);
+%! cest = condest (As);
+%! assert (cest, cexp, -1e-12);
+
 ## Test input validation
 %!error condest ()
 %!error condest (1,2,3,4,5,6,7)