[Octave-bug-tracker] [bug #54303] expm([]) results in error rather than []

[anonymous]

URL:
  <http://savannah.gnu.org/bugs/?54303>

                 Summary: expm([]) results in error rather than []
                 Project: GNU Octave
            Submitted by: None
            Submitted on: Fri 13 Jul 2018 07:51:10 AM UTC
                Category: Octave Function
                Severity: 3 - Normal
                Priority: 5 - Normal
              Item Group: Unexpected Error
                  Status: None
             Assigned to: None
         Originator Name: 
        Originator Email: address@hidden
             Open/Closed: Open
         Discussion Lock: Any
                 Release: 4.4.0
        Operating System: GNU/Linux

    _______________________________________________________

Details:

Functions like "sqrtm" and "logm" allow empty matrix as input.
But "expm" does not.
The expression "scalar^[]" also raises error.

octave:513>  logm ([])
ans = [](0x0)
octave:514>  sqrtm ([])
ans = [](0x0)
octave:515>  expm([])
warning: division by zero
warning: called from
    expm at line 96 column 11
 ** On entry to DGEBAL parameter number  4 had an illegal value
error: balance: unknown error in fortran subroutine
error: called from
    expm at line 101 column 11
octave:515>  3^[]
error: for x^A, A must be a square matrix.  Use .^ for elementwise power.



I'm not sure whether this is a bug.

If this is a bug, below is a possible fix for "expm.m".
(In this patch, only the first part is to fix this bug. The others are just
some small improvement.)

--- a/scripts/linear-algebra/expm.m
+++ b/scripts/linear-algebra/expm.m
@@ -85,17 +85,21 @@
   if (isscalar (A))
     r = exp (A);
     return;
+  elseif (isempty (A))
+    r = A;
+    return;
   elseif (strfind (typeinfo (A), "diagonal matrix"))
     r = diag (exp (diag (A)));
     return;
   endif
 
   n = rows (A);
+  id = eye (n);
   ## Trace reduction.
   A(A == -Inf) = -realmax;
-  trshift = trace (A) / length (A);
+  trshift = trace (A) / n;
   if (trshift > 0)
-    A -= trshift*eye (n);
+    A -= trshift * id;
   endif
   ## Balancing.
   [d, p, aa] = balance (A);
@@ -120,7 +124,6 @@
        1.9270852604185938e-9];
 
   a2 = aa^2;
-  id = eye (n);
   x = (((c(8) * a2 + c(6) * id) * a2 + c(4) * id) * a2 + c(2) * id) * a2 +
id;
   y = (((c(7) * a2 + c(5) * id) * a2 + c(3) * id) * a2 + c(1) * id) * aa;
 





    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?54303>

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