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
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
|
.TH DHGEQZ l "15 June 2000" "LAPACK version 3.0" ")"
.SH NAME
DHGEQZ - implement a single-/double-shift version of the QZ method for finding the generalized eigenvalues w(j)=(ALPHAR(j) + i*ALPHAI(j))/BETAR(j) of the equation det( A - w(i) B ) = 0 In addition, the pair A,B may be reduced to generalized Schur form
.SH SYNOPSIS
.TP 19
SUBROUTINE DHGEQZ(
JOB, COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB,
ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK,
LWORK, INFO )
.TP 19
.ti +4
CHARACTER
COMPQ, COMPZ, JOB
.TP 19
.ti +4
INTEGER
IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, LWORK, N
.TP 19
.ti +4
DOUBLE
PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
B( LDB, * ), BETA( * ), Q( LDQ, * ), WORK( * ),
Z( LDZ, * )
.SH PURPOSE
DHGEQZ implements a single-/double-shift version of the QZ method for finding the generalized eigenvalues w(j)=(ALPHAR(j) + i*ALPHAI(j))/BETAR(j) of the equation det( A - w(i) B ) = 0 In addition, the pair A,B may be reduced to generalized Schur form: B is upper triangular, and A is block upper triangular, where the
diagonal blocks are either 1-by-1 or 2-by-2, the 2-by-2 blocks having
complex generalized eigenvalues (see the description of the argument
JOB.)
.br
If JOB='S', then the pair (A,B) is simultaneously reduced to Schur
form by applying one orthogonal tranformation (usually called Q) on
the left and another (usually called Z) on the right. The 2-by-2
upper-triangular diagonal blocks of B corresponding to 2-by-2 blocks
of A will be reduced to positive diagonal matrices. (I.e.,
if A(j+1,j) is non-zero, then B(j+1,j)=B(j,j+1)=0 and B(j,j) and
B(j+1,j+1) will be positive.)
.br
If JOB='E', then at each iteration, the same transformations
are computed, but they are only applied to those parts of A and B
which are needed to compute ALPHAR, ALPHAI, and BETAR.
.br
If JOB='S' and COMPQ and COMPZ are 'V' or 'I', then the orthogonal
transformations used to reduce (A,B) are accumulated into the arrays
Q and Z s.t.:
.br
Q(in) A(in) Z(in)* = Q(out) A(out) Z(out)*
.br
Q(in) B(in) Z(in)* = Q(out) B(out) Z(out)*
.br
Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix
Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973),
pp. 241--256.
.br
.SH ARGUMENTS
.TP 8
JOB (input) CHARACTER*1
= 'E': compute only ALPHAR, ALPHAI, and BETA. A and B will
not necessarily be put into generalized Schur form.
= 'S': put A and B into generalized Schur form, as well
as computing ALPHAR, ALPHAI, and BETA.
.TP 8
COMPQ (input) CHARACTER*1
= 'N': do not modify Q.
.br
= 'V': multiply the array Q on the right by the transpose of
the orthogonal tranformation that is applied to the
left side of A and B to reduce them to Schur form.
= 'I': like COMPQ='V', except that Q will be initialized to
the identity first.
.TP 8
COMPZ (input) CHARACTER*1
= 'N': do not modify Z.
.br
= 'V': multiply the array Z on the right by the orthogonal
tranformation that is applied to the right side of
A and B to reduce them to Schur form.
= 'I': like COMPZ='V', except that Z will be initialized to
the identity first.
.TP 8
N (input) INTEGER
The order of the matrices A, B, Q, and Z. N >= 0.
.TP 8
ILO (input) INTEGER
IHI (input) INTEGER
It is assumed that A is already upper triangular in rows and
columns 1:ILO-1 and IHI+1:N.
1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
.TP 8
A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
On entry, the N-by-N upper Hessenberg matrix A. Elements
below the subdiagonal must be zero.
If JOB='S', then on exit A and B will have been
simultaneously reduced to generalized Schur form.
If JOB='E', then on exit A will have been destroyed.
The diagonal blocks will be correct, but the off-diagonal
portion will be meaningless.
.TP 8
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max( 1, N ).
.TP 8
B (input/output) DOUBLE PRECISION array, dimension (LDB, N)
On entry, the N-by-N upper triangular matrix B. Elements
below the diagonal must be zero. 2-by-2 blocks in B
corresponding to 2-by-2 blocks in A will be reduced to
positive diagonal form. (I.e., if A(j+1,j) is non-zero,
then B(j+1,j)=B(j,j+1)=0 and B(j,j) and B(j+1,j+1) will be
positive.)
If JOB='S', then on exit A and B will have been
simultaneously reduced to Schur form.
If JOB='E', then on exit B will have been destroyed.
Elements corresponding to diagonal blocks of A will be
correct, but the off-diagonal portion will be meaningless.
.TP 8
LDB (input) INTEGER
The leading dimension of the array B. LDB >= max( 1, N ).
.TP 8
ALPHAR (output) DOUBLE PRECISION array, dimension (N)
ALPHAR(1:N) will be set to real parts of the diagonal
elements of A that would result from reducing A and B to
Schur form and then further reducing them both to triangular
form using unitary transformations s.t. the diagonal of B
was non-negative real. Thus, if A(j,j) is in a 1-by-1 block
(i.e., A(j+1,j)=A(j,j+1)=0), then ALPHAR(j)=A(j,j).
Note that the (real or complex) values
(ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the
generalized eigenvalues of the matrix pencil A - wB.
.TP 8
ALPHAI (output) DOUBLE PRECISION array, dimension (N)
ALPHAI(1:N) will be set to imaginary parts of the diagonal
elements of A that would result from reducing A and B to
Schur form and then further reducing them both to triangular
form using unitary transformations s.t. the diagonal of B
was non-negative real. Thus, if A(j,j) is in a 1-by-1 block
(i.e., A(j+1,j)=A(j,j+1)=0), then ALPHAR(j)=0.
Note that the (real or complex) values
(ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the
generalized eigenvalues of the matrix pencil A - wB.
.TP 8
BETA (output) DOUBLE PRECISION array, dimension (N)
BETA(1:N) will be set to the (real) diagonal elements of B
that would result from reducing A and B to Schur form and
then further reducing them both to triangular form using
unitary transformations s.t. the diagonal of B was
non-negative real. Thus, if A(j,j) is in a 1-by-1 block
(i.e., A(j+1,j)=A(j,j+1)=0), then BETA(j)=B(j,j).
Note that the (real or complex) values
(ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the
generalized eigenvalues of the matrix pencil A - wB.
(Note that BETA(1:N) will always be non-negative, and no
BETAI is necessary.)
.TP 8
Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
If COMPQ='N', then Q will not be referenced.
If COMPQ='V' or 'I', then the transpose of the orthogonal
transformations which are applied to A and B on the left
will be applied to the array Q on the right.
.TP 8
LDQ (input) INTEGER
The leading dimension of the array Q. LDQ >= 1.
If COMPQ='V' or 'I', then LDQ >= N.
.TP 8
Z (input/output) DOUBLE PRECISION array, dimension (LDZ, N)
If COMPZ='N', then Z will not be referenced.
If COMPZ='V' or 'I', then the orthogonal transformations
which are applied to A and B on the right will be applied
to the array Z on the right.
.TP 8
LDZ (input) INTEGER
The leading dimension of the array Z. LDZ >= 1.
If COMPZ='V' or 'I', then LDZ >= N.
.TP 8
WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
.TP 8
LWORK (input) INTEGER
The dimension of the array WORK. LWORK >= max(1,N).
If LWORK = -1, then a workspace query is assumed; the routine
only calculates the optimal size of the WORK array, returns
this value as the first entry of the WORK array, and no error
message related to LWORK is issued by XERBLA.
.TP 8
INFO (output) INTEGER
= 0: successful exit
.br
< 0: if INFO = -i, the i-th argument had an illegal value
.br
= 1,...,N: the QZ iteration did not converge. (A,B) is not
in Schur form, but ALPHAR(i), ALPHAI(i), and
BETA(i), i=INFO+1,...,N should be correct.
= N+1,...,2*N: the shift calculation failed. (A,B) is not
in Schur form, but ALPHAR(i), ALPHAI(i), and
BETA(i), i=INFO-N+1,...,N should be correct.
> 2*N: various "impossible" errors.
.SH FURTHER DETAILS
Iteration counters:
.br
JITER -- counts iterations.
.br
IITER -- counts iterations run since ILAST was last
.br
changed. This is therefore reset only when a 1-by-1 or
2-by-2 block deflates off the bottom.
.br
|
