# HG changeset patch
# User Marco Caliari <marco.caliari@univr.it>
# Date 1331636195 -3600
# Node ID d2bffa78730e3862d5ebf538072b4a89f2cc6f52
# Parent  3d628878e1093c4eeec5598064718e53380c6396
Fix logm for complex matrix with real eigenvalues (bug #34893).

* crsf2csf, zrsf2csf: Fix off-by-one error.
* logm.m: Only truncate imaginary parts for real matrices. Add a test.
* schur.cc: Add a test for rsf2csf.x

diff --git a/libcruft/lapack-xtra/crsf2csf.f b/libcruft/lapack-xtra/crsf2csf.f
--- a/libcruft/lapack-xtra/crsf2csf.f
+++ b/libcruft/lapack-xtra/crsf2csf.f
@@ -25,6 +25,9 @@
        real c(n-1),s(n-1)
        real x,y,z
        integer j
+       do j = 1,n-1
+          c(j) = 1
+       end do
        j = 1
        do while (j < n)
 c apply previous rotations to rows
@@ -33,10 +36,9 @@
          y = t(j+1,j)
          if (y /= 0) then
 c 2x2 block, form Givens rotation [c, i*s; i*s, c]
-           x = t(j,j)
            z = t(j,j+1)
            c(j) = sqrt(z/(z-y))
-           s(j) = sign(sqrt(-y/(z-y)),z)
+           s(j) = sqrt(y/(y-z))
 c apply new rotation to t(j:j+1,j)
            call crcrot1(2,t(j,j),c(j),s(j))
 c apply all rotations to t(1:j+1,j+1)
@@ -45,10 +47,8 @@
            call crcrot2(j+1,t(1,j),t(1,j+1),c(j),s(j))
 c zero subdiagonal entry, skip next row
            t(j+1,j) = 0
-           c(j+1) = 1
            j = j + 2
          else
-           c(j) = 1
            j = j + 1
          end if
        end do
diff --git a/libcruft/lapack-xtra/zrsf2csf.f b/libcruft/lapack-xtra/zrsf2csf.f
--- a/libcruft/lapack-xtra/zrsf2csf.f
+++ b/libcruft/lapack-xtra/zrsf2csf.f
@@ -25,6 +25,9 @@
        double precision c(n-1),s(n-1)
        double precision x,y,z
        integer j
+       do j = 1,n-1
+          c(j) = 1
+       end do
        j = 1
        do while (j < n)
 c apply previous rotations to rows
@@ -33,10 +36,9 @@
          y = t(j+1,j)
          if (y /= 0) then
 c 2x2 block, form Givens rotation [c, i*s; i*s, c]
-           x = t(j,j)
            z = t(j,j+1)
            c(j) = sqrt(z/(z-y))
-           s(j) = sign(sqrt(-y/(z-y)),z)
+           s(j) = sqrt(y/(y-z))
 c apply new rotation to t(j:j+1,j)
            call zrcrot1(2,t(j,j),c(j),s(j))
 c apply all rotations to t(1:j+1,j+1)
@@ -45,10 +47,8 @@
            call zrcrot2(j+1,t(1,j),t(1,j+1),c(j),s(j))
 c zero subdiagonal entry, skip next row
            t(j+1,j) = 0
-           c(j+1) = 1
            j = j + 2
          else
-           c(j) = 1
            j = j + 1
          end if
        end do
diff --git a/scripts/linear-algebra/logm.m b/scripts/linear-algebra/logm.m
--- a/scripts/linear-algebra/logm.m
+++ b/scripts/linear-algebra/logm.m
@@ -104,7 +104,7 @@
   s = 2^k * u * s * u';
 
   ## Remove small complex values (O(eps)) which may have entered calculation
-  if (real_eig)
+  if (real_eig && isreal(A))
     s = real (s);
   endif
 
@@ -159,6 +159,7 @@
 %!assert(norm(logm([1 -1;0 1]) - [0 -1; 0 0]) < 1e-5);
 %!assert(norm(expm(logm([-1 2 ; 4 -1])) - [-1 2 ; 4 -1]) < 1e-5);
 %!assert(logm([1 -1 -1;0 1 -1; 0 0 1]), [0 -1 -1.5; 0 0 -1; 0 0 0], 1e-5);
+%!assert (logm (expm ([0 1i; -1i 0])), [0 1i; -1i 0], 10 * eps)
 
 %% Test input validation
 %!error logm ();
diff --git a/src/DLD-FUNCTIONS/schur.cc b/src/DLD-FUNCTIONS/schur.cc
--- a/src/DLD-FUNCTIONS/schur.cc
+++ b/src/DLD-FUNCTIONS/schur.cc
@@ -376,4 +376,9 @@
 %! assert (norm (tril (T, -1)), 0)
 %! assert (norm (U * U'), 1, 1e-14)
 
+%!test
+%! A = [0, 1;-1, 0];
+%! [u, t] = schur (A);
+%! [U, T] = rsf2csf (u,t);
+%! assert (U * T * U', A, 1e-14)
 */