Mercurial > hg > octave-lyh
annotate 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 |
| rev | line source |
|---|---|
|
14138
72c96de7a403
maint: update copyright notices for 2012
John W. Eaton <jwe@octave.org>
parents:
13141
diff
changeset
|
1 c Copyright (C) 2010-2012 VZLU Prague, a.s., Czech Republic |
| 10822 | 2 c |
| 3 c Author: Jaroslav Hajek <highegg@gmail.com> | |
| 4 c | |
| 5 c This file is part of Octave. | |
| 6 c | |
| 7 c Octave is free software; you can redistribute it and/or modify | |
| 8 c it under the terms of the GNU General Public License as published by | |
| 9 c the Free Software Foundation; either version 3 of the License, or | |
| 10 c (at your option) any later version. | |
| 11 c | |
| 12 c This program is distributed in the hope that it will be useful, | |
| 13 c but WITHOUT ANY WARRANTY; without even the implied warranty of | |
| 14 c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
| 15 c GNU General Public License for more details. | |
| 16 c | |
| 17 c You should have received a copy of the GNU General Public License | |
| 18 c along with this software; see the file COPYING. If not, see | |
| 19 c <http://www.gnu.org/licenses/>. | |
| 20 c | |
| 21 | |
| 22 subroutine crsf2csf(n,t,u,c,s) | |
| 23 integer n | |
| 24 complex t(n,n),u(n,n) | |
| 25 real c(n-1),s(n-1) | |
| 26 real x,y,z | |
| 27 integer j | |
|
14481
d2bffa78730e
Fix logm for complex matrix with real eigenvalues (bug #34893).
Marco Caliari <marco.caliari@univr.it>
parents:
14138
diff
changeset
|
28 do j = 1,n-1 |
|
d2bffa78730e
Fix logm for complex matrix with real eigenvalues (bug #34893).
Marco Caliari <marco.caliari@univr.it>
parents:
14138
diff
changeset
|
29 c(j) = 1 |
|
d2bffa78730e
Fix logm for complex matrix with real eigenvalues (bug #34893).
Marco Caliari <marco.caliari@univr.it>
parents:
14138
diff
changeset
|
30 end do |
| 10822 | 31 j = 1 |
| 32 do while (j < n) | |
| 33 c apply previous rotations to rows | |
| 34 call crcrot1(j,t(1,j),c,s) | |
| 35 | |
| 36 y = t(j+1,j) | |
| 37 if (y /= 0) then | |
|
13141
e81ddf9cacd5
maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
38 c 2x2 block, form Givens rotation [c, i*s; i*s, c] |
| 10822 | 39 z = t(j,j+1) |
|
10823
3d89d262f5d4
slight simplification in rsf2csf code
Jaroslav Hajek <highegg@gmail.com>
parents:
10822
diff
changeset
|
40 c(j) = sqrt(z/(z-y)) |
|
14481
d2bffa78730e
Fix logm for complex matrix with real eigenvalues (bug #34893).
Marco Caliari <marco.caliari@univr.it>
parents:
14138
diff
changeset
|
41 s(j) = sqrt(y/(y-z)) |
| 10822 | 42 c apply new rotation to t(j:j+1,j) |
| 43 call crcrot1(2,t(j,j),c(j),s(j)) | |
| 44 c apply all rotations to t(1:j+1,j+1) | |
| 45 call crcrot1(j+1,t(1,j+1),c,s) | |
| 46 c apply new rotation to columns j,j+1 | |
|
13141
e81ddf9cacd5
maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
47 call crcrot2(j+1,t(1,j),t(1,j+1),c(j),s(j)) |
| 10822 | 48 c zero subdiagonal entry, skip next row |
| 49 t(j+1,j) = 0 | |
| 50 j = j + 2 | |
| 51 else | |
| 52 j = j + 1 | |
| 53 end if | |
| 54 end do | |
| 55 | |
| 56 c apply rotations to last column if needed | |
| 57 if (j == n) then | |
| 58 call crcrot1(j,t(1,j),c,s) | |
| 59 end if | |
| 60 | |
| 61 c apply stored rotations to all columns of u | |
| 62 do j = 1,n-1 | |
| 63 if (c(j) /= 1) then | |
| 64 call crcrot2(n,u(1,j),u(1,j+1),c(j),s(j)) | |
| 65 end if | |
| 66 end do | |
| 67 | |
| 68 end subroutine | |
| 69 | |
| 70 subroutine crcrot1(n,x,c,s) | |
| 71 c apply rotations to a column from the left | |
| 72 integer n | |
| 73 complex x(n), t | |
| 74 real c(n-1),s(n-1) | |
| 75 integer i | |
| 76 do i = 1,n-1 | |
| 77 if (c(i) /= 1) then | |
| 78 t = x(i)*c(i) - x(i+1)*cmplx(0,s(i)) | |
| 79 x(i+1) = x(i+1)*c(i) - x(i)*cmplx(0,s(i)) | |
| 80 x(i) = t | |
| 81 endif | |
| 82 end do | |
| 83 end subroutine | |
| 84 | |
| 85 subroutine crcrot2(n,x,y,c,s) | |
| 86 c apply a single rotation from the right to a pair of columns | |
| 87 integer n | |
| 88 complex x(n),y(n),t | |
| 89 real c, s | |
| 90 integer i | |
| 91 do i = 1,n | |
| 92 t = x(i)*c + y(i)*cmplx(0,s) | |
| 93 y(i) = y(i)*c + x(i)*cmplx(0,s) | |
| 94 x(i) = t | |
| 95 end do | |
| 96 end subroutine | |
| 97 | |
| 98 | |
| 99 | |
| 100 | |
| 101 |