Mercurial > octave

changeset 31234:a1318deb4584

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080) sparse-xpow.cc: Return 0x0 matrix for empty inputs of all sizes. Return early for all-zeros matrix with positive power. Add FIXME for all-zeros matrix with negative power. Add BISTs.
author Arun Giridhar <arungiridhar@gmail.com>
date Tue, 20 Sep 2022 09:00:02 -0400
parents 95fa3182f00c
children b542b88ad3b6
files libinterp/corefcn/sparse-xpow.cc
diffstat 1 files changed, 43 insertions(+), 4 deletions(-) [+]
line wrap: on
line diff
--- a/libinterp/corefcn/sparse-xpow.cc	Mon Sep 19 13:14:16 2022 +0200
+++ b/libinterp/corefcn/sparse-xpow.cc	Tue Sep 20 09:00:02 2022 -0400
@@ -66,12 +66,13 @@
   octave_idx_type nr = a.rows ();
   octave_idx_type nc = a.cols ();
 
+  if (nr == 0 || nc == 0)
+    return SparseMatrix();
+
+  // If we are here, A is not empty ==> A needs to be square.
   if (nr != nc)
     error ("for A^b, A must be a square matrix.  Use .^ for elementwise power.");
 
-  if (nr == 0 && nc == 0)
-    return a;
-
   if (! xisint (b))
     error ("use full(a) ^ full(b)");
 
@@ -100,6 +101,9 @@
           double rcond = 0.0;
           MatrixType mattyp (a);
 
+          // FIXME: This causes an error if the input sparse matrix is all-zeros.
+          // That behavior is inconsistent with A ^ b when A is a full all-zeros
+          // matrix, which just returns Inf of the same size with a warning.
           atmp = a.inverse (mattyp, info, rcond, 1);
 
           if (info == -1)
@@ -108,6 +112,9 @@
       else
         atmp = a;
 
+      if (atmp.nnz() == 0)  // Fast return for all-zeros matrix
+        return atmp;
+
       SparseMatrix result (atmp);
 
       btmp--;
@@ -172,7 +179,11 @@
   octave_idx_type nr = a.rows ();
   octave_idx_type nc = a.cols ();
 
-  if (nr == 0 || nc == 0 || nr != nc)
+  if (nr == 0 || nc == 0)
+    return SparseMatrix();
+
+  // If we are here, A is not empty ==> A needs to be square.
+  if (nr != nc)
     error ("for A^b, A must be a square matrix.  Use .^ for elementwise power.");
 
   if (! xisint (b))
@@ -211,6 +222,9 @@
       else
         atmp = a;
 
+      if (atmp.nnz() == 0)  // Fast return for all-zeros matrix
+        return atmp;
+
       SparseComplexMatrix result (atmp);
 
       btmp--;
@@ -301,6 +315,30 @@
 /*
 %!assert (sparse (2) .^ [3, 4], sparse ([8, 16]))
 %!assert <47775> (sparse (2i) .^ [3, 4], sparse ([-0-8i, 16]))
+
+%!test <*63080>
+%! Z = sparse ([]);
+%! A = sparse (zeros (0, 2));
+%! B = sparse (zeros (2, 0));
+%! assert (Z ^  1, Z);
+%! assert (Z ^  0, Z);
+%! assert (Z ^ -1, Z);
+%! assert (A ^  1, Z);
+%! assert (A ^  0, Z);
+%! assert (A ^ -1, Z);
+%! assert (B ^  1, Z);
+%! assert (B ^  0, Z);
+%! assert (B ^ -1, Z);
+
+%!test <*63080>
+%! A = sparse (zeros (2, 2));
+%! assert (A ^  1, A);
+%! assert (A ^  0, sparse (eye (2, 2)));
+
+%!test <63080>
+%! A = sparse (zeros (2, 2));
+%! assert (A ^ -1, sparse (inf (2, 2)));
+
 */
 
 // -*- 1 -*-
@@ -803,4 +841,5 @@
   return result;
 }
 
+
 OCTAVE_NAMESPACE_END