octave-bug-tracker
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Octave-bug-tracker] [bug #56232] Octave crash when inverting an empty s


From: Rik
Subject: [Octave-bug-tracker] [bug #56232] Octave crash when inverting an empty sparse matrix.
Date: Sat, 18 May 2019 17:45:42 -0400 (EDT)
User-agent: Mozilla/5.0 (Windows NT 10.0; WOW64; Trident/7.0; rv:11.0) like Gecko

Update of bug #56232 (project octave):

                  Status:          Patch Reviewed => Fixed                  

    _______________________________________________________

Follow-up Comment #62:

Thanks to Marco for taking the lead in resolving this problem.  I made some
modifications for syntax or for performance, and checked it in here
(https://hg.savannah.gnu.org/hgweb/octave/rev/6e18f0ce268c).

Documenting my additions, I changed the BIST tests which were using some side
effects of our implementation to work, to use more ordinary syntax.

Instead of 


+%!test
+%!warning <matrix singular> A = inv (sparse ([1, 2;0 ,0]));
+%! assert (A, sparse ([Inf, Inf; 0, 0]));


I wrote


+%!test
+%! fail ("A = inv (sparse ([1, 2;0 ,0]))", "warning", "matrix singular");
+%! assert (A, sparse ([Inf, Inf; 0, 0]));


For the diagonal matrices, the implementation relied on creating a local
buffer that could hold up to nnz() elements.  The algorithm also needed to
cycle over the entire number of diagonal elements before arriving at a
decision about whether the matrix was singular.  The initial implementation is
below.


info = 0;
OCTAVE_LOCAL_BUFFER (octave_idx_type, nzelem, len);
octave_idx_type count_nz = 0;
for (octave_idx_type i = 0; i < len; i++)
  {
    if (elem (i, i) == 0.0)
      retval.elem (i, i) = octave::numeric_limits<double>::Inf ();
    else
      {
        retval.elem (i, i) = 1.0 / elem (i, i);
        nzelem[count_nz] = i;
        count_nz++;
      }
  }
if (count_nz == 0)
  {
    (*current_liboctave_error_handler)
      ("inverse of the null matrix not defined");
  }
else if (count_nz < len)
  {
    info = -1;
    for (octave_idx_type i = 0; i < count_nz; i++)
      retval.elem (nzelem[i], nzelem[i]) = octave::numeric_limits<double>::Inf
();
  }


Instead, I found that it was possible with two count variables to eliminate
the OCTAVE_LOCAL_BUFFER and to short-circuit out of the loop as soon as a zero
is detected.  For a very small speed-up, I also used xelem() to access
elements where I could be certain that I was not going to go outside the
bounds of the array.


info = 0;
octave_idx_type len = r;        // alias for readability
octave_idx_type z_count  = 0;   // zeros
octave_idx_type nz_count = 0;   // non-zeros
for (octave_idx_type i = 0; i < len; i++)
  {
    if (xelem (i, i) == 0.0)
      {
        z_count++;
        if (nz_count > 0)
          break;
      }
    else
      {
        nz_count++;
        if (z_count > 0)
          break;
        retval.elem (i, i) = 1.0 / xelem (i, i);
      }
  }
if (nz_count == 0)
  {
    (*current_liboctave_error_handler)
      ("inverse of the null matrix not defined");
  }
else if (z_count > 0)
  {
    info = -1;
    element_type *data = retval.fortran_vec ();
    std::fill (data, data + len, octave::numeric_limits<double>::Inf ());
  }


For Sparse and CompleSparse arrays, the original additional code was


+          if (rcond == 0.)
+            {
+              for (octave_idx_type i = 0; i < cols () + 1; i++)
+                ret.cidx (i) = cidx (i);
+
+              for (octave_idx_type i = 0; i < nnz (); i++)
+                {
+                  ret.data (i) = octave::numeric_limits<double>::Inf ();
+                  ret.ridx (i) = ridx (i);
+                }
+
+              return ret;
+            }


Instead of using hand-rolled for loops, I decided on functions from the C++
standard library.


+          if (rcond == 0.0)
+            {
+              // Return all Inf matrix with sparsity pattern of input.
+              octave_idx_type nz = nnz ();
+              ret = SparseComplexMatrix (rows (), cols (), nz);
+              std::fill (ret.xdata (), ret.xdata () + nz,
+                         octave::numeric_limits<double>::Inf ());
+              std::copy_n (ridx (), nz, ret.xridx ());
+              std::copy_n (cidx (), cols () + 1, ret.xcidx ());
+
+              return ret;
+            }




    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?56232>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/




reply via email to

[Prev in Thread] Current Thread [Next in Thread]