annotate liboctave/external/lapack-xtra/zrsf2csf.f @ 30564:796f54d4ddbf stable

update Octave Project Developers copyright for the new year In files that have the "Octave Project Developers" copyright notice, update for 2021. In all .txi and .texi files except gpl.txi and gpl.texi in the doc/liboctave and doc/interpreter directories, change the copyright to "Octave Project Developers", the same as used for other source files. Update copyright notices for 2022 (not done since 2019). For gpl.txi and gpl.texi, change the copyright notice to be "Free Software Foundation, Inc." and leave the date at 2007 only because this file only contains the text of the GPL, not anything created by the Octave Project Developers. Add Paul Thomas to contributors.in.
author John W. Eaton <jwe@octave.org>
date Tue, 28 Dec 2021 18:22:40 -0500
parents 0a5b15007766
children 597f3ee61a48
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
30564
796f54d4ddbf update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents: 29358
diff changeset
1 c Copyright (C) 2010-2022 The Octave Project Developers
27918
b442ec6dda5c use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents: 26376
diff changeset
2 c
27923
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
3 c See the file COPYRIGHT.md in the top-level directory of this
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
4 c distribution or <https://octave.org/copyright/>.
10822
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
5 c
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
6 c This file is part of Octave.
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
7 c
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23434
diff changeset
8 c Octave is free software: you can redistribute it and/or modify it
22802
0dcff7695e26 maint: Update more Copyright statements to use standard form.
Rik <rik@octave.org>
parents: 22323
diff changeset
9 c under the terms of the GNU General Public License as published by
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23434
diff changeset
10 c the Free Software Foundation, either version 3 of the License, or
10822
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
11 c (at your option) any later version.
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
12 c
22802
0dcff7695e26 maint: Update more Copyright statements to use standard form.
Rik <rik@octave.org>
parents: 22323
diff changeset
13 c Octave is distributed in the hope that it will be useful, but
0dcff7695e26 maint: Update more Copyright statements to use standard form.
Rik <rik@octave.org>
parents: 22323
diff changeset
14 c WITHOUT ANY WARRANTY; without even the implied warranty of
10822
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
15 c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
16 c GNU General Public License for more details.
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
17 c
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
18 c You should have received a copy of the GNU General Public License
22802
0dcff7695e26 maint: Update more Copyright statements to use standard form.
Rik <rik@octave.org>
parents: 22323
diff changeset
19 c along with Octave; see the file COPYING. If not, see
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23434
diff changeset
20 c <https://www.gnu.org/licenses/>.
10822
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
21 c
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
22
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
23 subroutine zrsf2csf(n,t,u,c,s)
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
24 integer n
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
25 double complex t(n,n),u(n,n)
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
26 double precision c(n-1),s(n-1)
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
27 double precision x,y,z
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
28 integer j
14460
6c3441f3146b Fix logm for complex matrix with real eigenvalues (bug #34893).
Marco Caliari <marco.caliari@univr.it>
parents: 14138
diff changeset
29 do j = 1,n-1
6c3441f3146b Fix logm for complex matrix with real eigenvalues (bug #34893).
Marco Caliari <marco.caliari@univr.it>
parents: 14138
diff changeset
30 c(j) = 1
6c3441f3146b Fix logm for complex matrix with real eigenvalues (bug #34893).
Marco Caliari <marco.caliari@univr.it>
parents: 14138
diff changeset
31 end do
10822
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
32 j = 1
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
33 do while (j < n)
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
34 c apply previous rotations to rows
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
35 call zrcrot1(j,t(1,j),c,s)
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
36
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
37 y = t(j+1,j)
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
38 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
39 c 2x2 block, form Givens rotation [c, i*s; i*s, c]
10822
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
40 z = t(j,j+1)
10823
3d89d262f5d4 slight simplification in rsf2csf code
Jaroslav Hajek <highegg@gmail.com>
parents: 10822
diff changeset
41 c(j) = sqrt(z/(z-y))
14460
6c3441f3146b Fix logm for complex matrix with real eigenvalues (bug #34893).
Marco Caliari <marco.caliari@univr.it>
parents: 14138
diff changeset
42 s(j) = sqrt(y/(y-z))
10822
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
43 c apply new rotation to t(j:j+1,j)
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
44 call zrcrot1(2,t(j,j),c(j),s(j))
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
45 c apply all rotations to t(1:j+1,j+1)
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
46 call zrcrot1(j+1,t(1,j+1),c,s)
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
47 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
48 call zrcrot2(j+1,t(1,j),t(1,j+1),c(j),s(j))
10822
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
49 c zero subdiagonal entry, skip next row
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
50 t(j+1,j) = 0
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
51 j = j + 2
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
52 else
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
53 j = j + 1
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
54 end if
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
55 end do
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
56
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
57 c apply rotations to last column if needed
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
58 if (j == n) then
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
59 call zrcrot1(j,t(1,j),c,s)
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
60 end if
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
61
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
62 c apply stored rotations to all columns of u
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
63 do j = 1,n-1
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
64 if (c(j) /= 1) then
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
65 call zrcrot2(n,u(1,j),u(1,j+1),c(j),s(j))
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
66 end if
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
67 end do
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
68
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
69 end subroutine
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
70
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
71 subroutine zrcrot1(n,x,c,s)
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
72 c apply rotations to a column from the left
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
73 integer n
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
74 double complex x(n), t
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
75 double precision c(n-1),s(n-1)
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
76 integer i
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
77 do i = 1,n-1
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
78 if (c(i) /= 1) then
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
79 t = x(i)*c(i) - x(i+1)*dcmplx(0,s(i))
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
80 x(i+1) = x(i+1)*c(i) - x(i)*dcmplx(0,s(i))
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
81 x(i) = t
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
82 endif
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
83 end do
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
84 end subroutine
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
85
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
86 subroutine zrcrot2(n,x,y,c,s)
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
87 c apply a single rotation from the right to a pair of columns
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
88 integer n
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
89 double complex x(n),y(n),t
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
90 double precision c, s
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
91 integer i
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
92 do i = 1,n
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
93 t = x(i)*c + y(i)*dcmplx(0,s)
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
94 y(i) = y(i)*c + x(i)*dcmplx(0,s)
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
95 x(i) = t
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
96 end do
23d2378512a0 implement rsf2csf
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
97 end subroutine