GNU Octave  3.8.0 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
zrsf2csf.f
Go to the documentation of this file.
1 c Copyright (C) 2010-2013 VZLU Prague, a.s., Czech Republic
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
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
20 c
21
22  subroutine zrsf2csf(n,t,u,c,s)
23  integer n
24  double complex t(n,n),u(n,n)
25  double precision c(n-1),s(n-1)
26  double precision x,y,z
27  integer j
28  do j = 1,n-1
29  c(j) = 1
30  end do
31  j = 1
32  do while (j < n)
33 c apply previous rotations to rows
34  call zrcrot1(j,t(1,j),c,s)
35
36  y = t(j+1,j)
37  if (y /= 0) then
38 c 2x2 block, form Givens rotation [c, i*s; i*s, c]
39  z = t(j,j+1)
40  c(j) = sqrt(z/(z-y))
41  s(j) = sqrt(y/(y-z))
42 c apply new rotation to t(j:j+1,j)
43  call zrcrot1(2,t(j,j),c(j),s(j))
44 c apply all rotations to t(1:j+1,j+1)
45  call zrcrot1(j+1,t(1,j+1),c,s)
46 c apply new rotation to columns j,j+1
47  call zrcrot2(j+1,t(1,j),t(1,j+1),c(j),s(j))
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 zrcrot1(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 zrcrot2(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 zrcrot1(n,x,c,s)
71 c apply rotations to a column from the left
72  integer n
73  double complex x(n), t
74  double precision 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)*dcmplx(0,s(i))
79  x(i+1) = x(i+1)*c(i) - x(i)*dcmplx(0,s(i))
80  x(i) = t
81  endif
82  end do
83  end subroutine
84
85  subroutine zrcrot2(n,x,y,c,s)
86 c apply a single rotation from the right to a pair of columns
87  integer n
88  double complex x(n),y(n),t
89  double precision c, s
90  integer i
91  do i = 1,n
92  t = x(i)*c + y(i)*dcmplx(0,s)
93  y(i) = y(i)*c + x(i)*dcmplx(0,s)
94  x(i) = t
95  end do
96  end subroutine
97
98
99
100
101