Mercurial > octave

changeset 29777:90d7137b7dc6

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
sqrtm: Fix for input with eigenvalues 0 (bug #60797). * libinterp/corefcn/sqrtm.cc (sqrtm_utri_inplace): Do not divide elements that are 0 (to avoid NaN). Add test.
author Markus Mützel <markus.muetzel@gmx.de>
date Fri, 18 Jun 2021 16:09:10 +0200
parents 0916ffc997e6
children c44e54bb018b
files libinterp/corefcn/sqrtm.cc
diffstat 1 files changed, 7 insertions(+), 2 deletions(-) [+]
line wrap: on
line diff
--- a/libinterp/corefcn/sqrtm.cc	Fri Jun 18 10:11:06 2021 +0200
+++ b/libinterp/corefcn/sqrtm.cc	Fri Jun 18 16:09:10 2021 +0200
@@ -54,7 +54,9 @@
   //   for j = 1:n
   //     T(j,j) = sqrt (T(j,j));
   //     for i = j-1:-1:1
-  //       T(i,j) /= (T(i,i) + T(j,j));
+  //       if T(i,j) != 0
+  //         T(i,j) /= (T(i,i) + T(j,j));
+  //       endif
   //       k = 1:i-1;
   //       T(k,j) -= T(k,i) * T(i,j);
   //     endfor
@@ -76,7 +78,9 @@
       for (octave_idx_type i = j-1; i >= 0; i--)
         {
           const element_type *coli = Tp + n*i;
-          const element_type colji = colj[i] /= (coli[i] + colj[j]);
+          if (colj[i] != zero)
+            colj[i] /= (coli[i] + colj[j]);
+          const element_type colji = colj[i];
           for (octave_idx_type k = 0; k < i; k++)
             colj[k] -= coli[k] * colji;
         }
@@ -256,6 +260,7 @@
 
 /*
 %!assert (sqrtm (2*ones (2)), ones (2), 3*eps)
+%!assert (real (sqrtm (ones (4))), 0.5*ones (4), 4*eps)
 
 ## The following two tests are from the reference in the docstring above.
 %!test