Mercurial > hg > octave-lyh
changeset 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 | 3d628878e109 |
| children | 51fd0cf227e4 |
| files | libcruft/lapack-xtra/crsf2csf.f libcruft/lapack-xtra/zrsf2csf.f scripts/linear-algebra/logm.m src/DLD-FUNCTIONS/schur.cc |
| diffstat | 4 files changed, 15 insertions(+), 9 deletions(-) [+] |
line wrap: on
line diff
--- 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
--- 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
--- 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 ();