1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
|
SUBROUTINE ZLANV2( A, B, C, D, RT1, RT2, CS, SN )
*
* -- ScaLAPACK routine (version 1.7) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* May 28, 1999
*
* .. Scalar Arguments ..
DOUBLE PRECISION CS
COMPLEX*16 A, B, C, D, RT1, RT2, SN
* ..
*
* Purpose
* =======
*
* ZLANV2 computes the Schur factorization of a complex 2-by-2
* nonhermitian matrix in standard form:
*
* [ A B ] = [ CS -SN ] [ AA BB ] [ CS SN ]
* [ C D ] [ SN CS ] [ 0 DD ] [-SN CS ]
*
* Arguments
* =========
*
* A (input/output) COMPLEX*16
* B (input/output) COMPLEX*16
* C (input/output) COMPLEX*16
* D (input/output) COMPLEX*16
* On entry, the elements of the input matrix.
* On exit, they are overwritten by the elements of the
* standardised Schur form.
*
* RT1 (output) COMPLEX*16
* RT2 (output) COMPLEX*16
* The two eigenvalues.
*
* CS (output) DOUBLE PRECISION
* SN (output) COMPLEX*16
* Parameters of the rotation matrix.
*
* Further Details
* ===============
*
* Implemented by Mark R. Fahey, May 28, 1999
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION RZERO, HALF, RONE
PARAMETER ( RZERO = 0.0D+0, HALF = 0.5D+0,
$ RONE = 1.0D+0 )
COMPLEX*16 ZERO, ONE
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ),
$ ONE = ( 1.0D+0, 0.0D+0 ) )
* ..
* .. Local Scalars ..
COMPLEX*16 AA, BB, DD, T, TEMP, TEMP2, U, X, Y
* ..
* .. External Functions ..
COMPLEX*16 ZLADIV
EXTERNAL ZLADIV
* ..
* .. External Subroutines ..
EXTERNAL ZLARTG
* ..
* .. Intrinsic Functions ..
INTRINSIC DBLE, DCMPLX, DCONJG, DIMAG, SQRT
* ..
* .. Executable Statements ..
*
* Initialize CS and SN
*
CS = RONE
SN = ZERO
*
IF( C.EQ.ZERO ) THEN
GO TO 10
*
ELSE IF( B.EQ.ZERO ) THEN
*
* Swap rows and columns
*
CS = RZERO
SN = ONE
TEMP = D
D = A
A = TEMP
B = -C
C = ZERO
GO TO 10
ELSE IF( ( A-D ).EQ.ZERO ) THEN
TEMP = SQRT( B*C )
A = A + TEMP
D = D - TEMP
IF( ( B+C ).EQ.ZERO ) THEN
CS = SQRT( HALF )
SN = DCMPLX( RZERO, RONE )*CS
ELSE
TEMP = SQRT( B+C )
TEMP2 = ZLADIV( SQRT( B ), TEMP )
CS = DBLE( TEMP2 )
SN = ZLADIV( SQRT( C ), TEMP )
END IF
B = B - C
C = ZERO
GO TO 10
ELSE
*
* Compute eigenvalue closest to D
*
T = D
U = B*C
X = HALF*( A-T )
Y = SQRT( X*X+U )
IF( DBLE( X )*DBLE( Y )+DIMAG( X )*DIMAG( Y ).LT.RZERO )
$ Y = -Y
T = T - ZLADIV( U, ( X+Y ) )
*
* Do one QR step with exact shift T - resulting 2 x 2 in
* triangular form.
*
CALL ZLARTG( A-T, C, CS, SN, AA )
*
D = D - T
BB = CS*B + SN*D
DD = -DCONJG( SN )*B + CS*D
*
A = AA*CS + BB*DCONJG( SN ) + T
B = -AA*SN + BB*CS
C = ZERO
D = T
*
END IF
*
10 CONTINUE
*
* Store eigenvalues in RT1 and RT2.
*
RT1 = A
RT2 = D
RETURN
*
* End of ZLANV2
*
END
|
