Mercurial > octave-nkf
diff libcruft/lapack-xtra/crsf2csf.f @ 14481:d2bffa78730e stable
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
| author | Marco Caliari <marco.caliari@univr.it> |
|---|---|
| date | Tue, 13 Mar 2012 11:56:35 +0100 |
| parents | 72c96de7a403 |
| children |
line wrap: on
line diff
--- a/libcruft/lapack-xtra/crsf2csf.f Sat Mar 17 23:24:10 2012 -0400 +++ b/libcruft/lapack-xtra/crsf2csf.f Tue Mar 13 11:56:35 2012 +0100 @@ -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