Mercurial > octave

changeset 20707:cd1bd06974d8

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Preserve diagonal matrix property with linear index assignment (bug #36932). * ov-base-diag.cc (subsasgn): Add new if test to check for a single linear index. Use sub2ind to convert index to subscripts. If subscripts i,j are equal then assignment is to a diagonal and diagonal matrix is preserved. * test/diag-perm.tst: Add new tests to verify preservation of diagonal matrix property under subscript and index assignment.
author Rik <rik@octave.org>
date Mon, 16 Nov 2015 21:27:40 -0800
parents fec7cc73507b
children 453fca9ae397
files libinterp/octave-value/ov-base-diag.cc test/diag-perm.tst
diffstat 2 files changed, 57 insertions(+), 5 deletions(-) [+]
line wrap: on
line diff
--- a/libinterp/octave-value/ov-base-diag.cc	Mon Nov 02 22:42:28 2015 +0100
+++ b/libinterp/octave-value/ov-base-diag.cc	Mon Nov 16 21:27:40 2015 -0800
@@ -158,11 +158,46 @@
         if (type.length () == 1)
           {
             octave_value_list jdx = idx.front ();
-            // Check for a simple element assignment. That means, if D is a
-            // diagonal matrix, 'D(i,i) = x' will not destroy its diagonality
-            // (provided i is a valid index).
-            if (jdx.length () == 2
-                && jdx(0).is_scalar_type () && jdx(1).is_scalar_type ())
+
+            // FIXME: Mostly repeated code for cases 1 and 2 could be
+            //        consolidated for DRY (Don't Repeat Yourself).
+            // Check for assignments to diagonal elements which should not
+            // destroy the diagonal property of the matrix.
+            // If D is a diagonal matrix then the assignment can be
+            // 1) linear, D(i) = x, where ind2sub results in case #2 below
+            // 2) subscript D(i,i) = x, where both indices are equal.
+            if (jdx.length () == 1 && jdx(0).is_scalar_type ())
+              {
+                typename DMT::element_type val;
+                int k = 0;
+                try
+                  {
+                    idx_vector ind = jdx(0).index_vector ();
+                    k = 1;
+                    dim_vector dv (matrix.rows (), matrix.cols ());
+                    Array<idx_vector> ivec = ind2sub (dv, ind);
+                    idx_vector i0 = ivec(0);
+                    idx_vector i1 = ivec(1);
+
+                    if (i0(0) == i1(0)
+                        && chk_valid_scalar (rhs, val))
+                      {
+                        matrix.dgelem (i0(0)) = val;
+                        retval = this;
+                        this->count++;
+                        // invalidate cache
+                        dense_cache = octave_value ();
+                      }
+                  }
+                catch (index_exception& e)
+                  {
+                    // Rethrow to allow more info to be reported later.
+                    e.set_pos_if_unset (2, k+1);
+                    throw;
+                  }
+              }
+            else if (jdx.length () == 2
+                     && jdx(0).is_scalar_type () && jdx(1).is_scalar_type ())
               {
                 typename DMT::element_type val;
                 int k = 0;
--- a/test/diag-perm.tst	Mon Nov 02 22:42:28 2015 +0100
+++ b/test/diag-perm.tst	Mon Nov 16 21:27:40 2015 -0800
@@ -270,3 +270,20 @@
 %! assert (inv (x), diag ([1 1/2 1/3]));
 %! x = diag (0:2);
 %! assert (inv (x), diag ([Inf 1 1/2]));
+
+## assignment to diagonal elements preserves diagonal structure (bug #36932)
+%!test
+%! x = diag (1:3);
+%! x(1,1) = -1;
+%! assert (typeinfo (x), "diagonal matrix");
+%! x(3,3) = -1;
+%! assert (typeinfo (x), "diagonal matrix");
+
+%!test
+%! x = diag (1:3);
+%! x(1) = -1;
+%! assert (typeinfo (x), "diagonal matrix");
+%! x(9) = -1;
+%! assert (typeinfo (x), "diagonal matrix");
+
+