GNU Octave  6.2.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-2021 The Octave Project Developers
2 c
3 c See the file COPYRIGHT.md in the top-level directory of this
4 c distribution or <https://octave.org/copyright/>.
5 c
6 c This file is part of Octave.
7 c
8 c Octave is free software: you can redistribute it and/or modify it
9 c under the terms of the GNU General Public License as published by
10 c the Free Software Foundation, either version 3 of the License, or
11 c (at your option) any later version.
12 c
13 c Octave is distributed in the hope that it will be useful, but
14 c WITHOUT ANY WARRANTY; without even the implied warranty of
15 c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 c GNU General Public License for more details.
17 c
18 c You should have received a copy of the GNU General Public License
19 c along with Octave; see the file COPYING. If not, see
20 c <https://www.gnu.org/licenses/>.
21 c
22 
23  subroutine zrsf2csf(n,t,u,c,s)
24  integer n
25  double complex t(n,n),u(n,n)
26  double precision c(n-1),s(n-1)
27  double precision x,y,z
28  integer j
29  do j = 1,n-1
30  c(j) = 1
31  end do
32  j = 1
33  do while (j < n)
34 c apply previous rotations to rows
35  call zrcrot1(j,t(1,j),c,s)
36 
37  y = t(j+1,j)
38  if (y /= 0) then
39 c 2x2 block, form Givens rotation [c, i*s; i*s, c]
40  z = t(j,j+1)
41  c(j) = sqrt(z/(z-y))
42  s(j) = sqrt(y/(y-z))
43 c apply new rotation to t(j:j+1,j)
44  call zrcrot1(2,t(j,j),c(j),s(j))
45 c apply all rotations to t(1:j+1,j+1)
46  call zrcrot1(j+1,t(1,j+1),c,s)
47 c apply new rotation to columns j,j+1
48  call zrcrot2(j+1,t(1,j),t(1,j+1),c(j),s(j))
49 c zero subdiagonal entry, skip next row
50  t(j+1,j) = 0
51  j = j + 2
52  else
53  j = j + 1
54  end if
55  end do
56 
57 c apply rotations to last column if needed
58  if (j == n) then
59  call zrcrot1(j,t(1,j),c,s)
60  end if
61 
62 c apply stored rotations to all columns of u
63  do j = 1,n-1
64  if (c(j) /= 1) then
65  call zrcrot2(n,u(1,j),u(1,j+1),c(j),s(j))
66  end if
67  end do
68 
69  end subroutine
70 
71  subroutine zrcrot1(n,x,c,s)
72 c apply rotations to a column from the left
73  integer n
74  double complex x(n), t
75  double precision c(n-1),s(n-1)
76  integer i
77  do i = 1,n-1
78  if (c(i) /= 1) then
79  t = x(i)*c(i) - x(i+1)*dcmplx(0,s(i))
80  x(i+1) = x(i+1)*c(i) - x(i)*dcmplx(0,s(i))
81  x(i) = t
82  endif
83  end do
84  end subroutine
85 
86  subroutine zrcrot2(n,x,y,c,s)
87 c apply a single rotation from the right to a pair of columns
88  integer n
89  double complex x(n),y(n),t
90  double precision c, s
91  integer i
92  do i = 1,n
93  t = x(i)*c + y(i)*dcmplx(0,s)
94  y(i) = y(i)*c + x(i)*dcmplx(0,s)
95  x(i) = t
96  end do
97  end subroutine
subroutine zrcrot2(n, x, y, c, s)
Definition: zrsf2csf.f:87
subroutine zrcrot1(n, x, c, s)
Definition: zrsf2csf.f:72
subroutine zrsf2csf(n, t, u, c, s)
Definition: zrsf2csf.f:24