# HG changeset patch
# User Rik <rik@octave.org>
# Date 1583957498 25200
# Node ID 4609d001daeefb6939b923a78806b332c91091f5
# Parent  026bff6a54d7ed77dd04baa14d62692585086b4f
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.

diff -r 026bff6a54d7 -r 4609d001daee scripts/linear-algebra/condest.m
--- 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)