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
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
|
<HTML>
<HEAD><TITLE>MB04HD - SLICOT Library Routine Documentation</TITLE>
</HEAD>
<BODY>
<H2><A Name="MB04HD">MB04HD</A></H2>
<H3>
Reducing a special real block (anti-)diagonal skew-Hamiltonian/Hamiltonian pencil to generalized Schur form
</H3>
<A HREF ="#Specification"><B>[Specification]</B></A>
<A HREF ="#Arguments"><B>[Arguments]</B></A>
<A HREF ="#Method"><B>[Method]</B></A>
<A HREF ="#References"><B>[References]</B></A>
<A HREF ="#Comments"><B>[Comments]</B></A>
<A HREF ="#Example"><B>[Example]</B></A>
<P>
<B><FONT SIZE="+1">Purpose</FONT></B>
<PRE>
To compute the transformed matrices A and B, using orthogonal
matrices Q1 and Q2 for a real N-by-N regular pencil
( A11 0 ) ( 0 B12 )
aA - bB = a ( ) - b ( ), (1)
( 0 A22 ) ( B21 0 )
where A11, A22 and B12 are upper triangular, B21 is upper
quasi-triangular and the generalized matrix product
-1 -1
A11 B12 A22 B21 is in periodic Schur form, such that the
matrix Q2' A Q1 is upper triangular, Q2' B Q1 is upper
quasi-triangular and the transformed pencil
a(Q2' A Q1) - b(Q2' B Q1) is in generalized Schur form. The
notation M' denotes the transpose of the matrix M.
</PRE>
<A name="Specification"><B><FONT SIZE="+1">Specification</FONT></B></A>
<PRE>
SUBROUTINE MB04HD( COMPQ1, COMPQ2, N, A, LDA, B, LDB, Q1, LDQ1,
$ Q2, LDQ2, IWORK, LIWORK, DWORK, LDWORK, BWORK,
$ INFO )
C .. Scalar Arguments ..
CHARACTER COMPQ1, COMPQ2
INTEGER INFO, LDA, LDB, LDQ1, LDQ2, LDWORK, LIWORK, N
C .. Array Arguments ..
LOGICAL BWORK( * )
INTEGER IWORK( * )
DOUBLE PRECISION A( LDA, * ), B( LDB, * ), DWORK( * ),
$ Q1( LDQ1, * ), Q2( LDQ2, * )
</PRE>
<A name="Arguments"><B><FONT SIZE="+1">Arguments</FONT></B></A>
<P>
<B>Mode Parameters</B>
<PRE>
COMPQ1 CHARACTER*1
Specifies whether to compute the orthogonal transformation
matrix Q1, as follows:
= 'N': Q1 is not computed;
= 'I': the array Q1 is initialized internally to the unit
matrix, and the orthogonal matrix Q1 is returned;
= 'U': the array Q1 contains an orthogonal matrix Q01 on
entry, and the matrix Q01*Q1 is returned, where Q1
is the product of the orthogonal transformations
that are applied on the right to the pencil
aA - bB in (1).
COMPQ2 CHARACTER*1
Specifies whether to compute the orthogonal transformation
matrix Q2, as follows:
= 'N': Q2 is not computed;
= 'I': the array Q2 is initialized internally to the unit
matrix, and the orthogonal matrix Q2 is returned;
= 'U': the array Q2 contains an orthogonal matrix Q02 on
entry, and the matrix Q02*Q2 is returned, where Q2
is the product of the orthogonal transformations
that are applied on the left to the pencil aA - bB
in (1).
</PRE>
<B>Input/Output Parameters</B>
<PRE>
N (input) INTEGER
Order of the pencil aA - bB, N >= 0, even.
A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
On entry, the leading N-by-N block diagonal part of this
array must contain the matrix A in (1). The off-diagonal
blocks need not be set to zero.
On exit, the leading N-by-N part of this array contains
the transformed upper triangular matrix.
LDA INTEGER
The leading dimension of the array A. LDA >= MAX(1, N).
B (input/output) DOUBLE PRECISION array, dimension (LDB, N)
On entry, the leading N-by-N block anti-diagonal part of
this array must contain the matrix B in (1). The diagonal
blocks need not be set to zero.
On exit, the leading N-by-N part of this array contains
the transformed upper quasi-triangular matrix.
LDB INTEGER
The leading dimension of the array B. LDB >= MAX(1, N).
Q1 (input/output) DOUBLE PRECISION array, dimension (LDQ1, N)
On entry, if COMPQ1 = 'U', then the leading N-by-N part of
this array must contain a given matrix Q01, and on exit,
the leading N-by-N part of this array contains the product
of the input matrix Q01 and the transformation matrix Q1
used to transform the matrices A and B.
On exit, if COMPQ1 = 'I', then the leading N-by-N part of
this array contains the orthogonal transformation matrix
Q1.
If COMPQ1 = 'N' this array is not referenced.
LDQ1 INTEGER
The leading dimension of the array Q1.
LDQ1 >= 1, if COMPQ1 = 'N';
LDQ1 >= MAX(1, N), if COMPQ1 = 'I' or COMPQ1 = 'U'.
Q2 (input/output) DOUBLE PRECISION array, dimension (LDQ2, N)
On entry, if COMPQ2 = 'U', then the leading N-by-N part of
this array must contain a given matrix Q02, and on exit,
the leading N-by-N part of this array contains the product
of the input matrix Q02 and the transformation matrix Q2
used to transform the matrices A and B.
On exit, if COMPQ2 = 'I', then the leading N-by-N part of
this array contains the orthogonal transformation matrix
Q2.
If COMPQ2 = 'N' this array is not referenced.
LDQ2 INTEGER
The leading dimension of the array Q2.
LDQ2 >= 1, if COMPQ2 = 'N';
LDQ2 >= MAX(1, N), if COMPQ2 = 'I' or COMPQ2 = 'U'.
</PRE>
<B>Workspace</B>
<PRE>
IWORK INTEGER array, dimension (LIWORK)
LIWORK INTEGER
The dimension of the array IWORK.
LIWORK >= MAX( N/2+1, 32 ).
DWORK DOUBLE PRECISION array, dimension (LDWORK)
On exit, if INFO = 0, DWORK(1) returns the optimal LDWORK.
On exit, if INFO = -16, DWORK(1) returns the minimum value
of LDWORK.
LDWORK INTEGER
The dimension of the array DWORK.
LDWORK >= 2*N*N + MAX( N/2 + 168, 272 ).
For good performance LDWORK should be generally larger.
If LDWORK = -1 a workspace query is assumed; the
routine only calculates the optimal size of the DWORK
array, returns this value as the first entry of the DWORK
array, and no error message is issued by XERBLA.
BWORK LOGICAL array, dimension (N/2)
</PRE>
<B>Error Indicator</B>
<PRE>
INFO INTEGER
= 0: succesful exit;
< 0: if INFO = -i, the i-th argument had an illegal value;
= 1: the periodic QZ algorithm failed to reorder the
eigenvalues (the problem is very ill-conditioned) in
the SLICOT Library routine MB03KD;
= 2: the standard QZ algorithm failed in the LAPACK
routines DGGES or DHGEQZ, called by the SLICOT
routines MB03DD or MB03FD;
= 3: the eigenvalue reordering failed in the LAPACK
routine DTGEX2, called by the SLICOT routine MB03FD;
= 4: the standard QZ algorithm failed to reorder the
eigenvalues in the LAPACK routine DTGSEN, called by
the SLICOT routine MB03DD.
</PRE>
<A name="Method"><B><FONT SIZE="+1">Method</FONT></B></A>
<PRE>
First, the periodic QZ algorithm (see also [2] and [3]) is applied
-1 -1
to the formal matrix product A11 B12 A22 B21 to reorder the
eigenvalues, i.e., orthogonal matrices V1, V2, V3 and V4 are
computed such that V2' A11 V1, V2' B12 V3, V4' A22 V3 and
V4' B21 V1 keep the triangular form, but they can be partitioned
into 2-by-2 block forms and the last diagonal blocks correspond to
all nonpositive real eigenvalues of the formal product, and the
first diagonal blocks correspond to the remaining eigenvalues.
Second, Q1 = diag(V1, V3), Q2 = diag(V2, V4) and
( AA11 AA12 0 0 )
( )
( 0 AA22 0 0 )
A := Q2' A Q1 =: ( ),
( 0 0 AA33 AA34 )
( )
( 0 0 0 AA44 )
( 0 0 BB13 BB14 )
( )
( 0 0 0 BB24 )
B := Q2' B Q1 =: ( ),
( BB31 BB32 0 0 )
( )
( 0 BB42 0 0 )
-1 -1
are set, such that AA22 BB24 AA44 BB42 has only nonpositive
real eigenvalues.
Third, the permutation matrix
( I 0 0 0 )
( )
( 0 0 I 0 )
P = ( ),
( 0 I 0 0 )
( )
( 0 0 0 I )
where I denotes the identity matrix of appropriate size, is used
to transform aA - bB to block upper triangular form
( AA11 0 | AA12 0 )
( | )
( 0 AA33 | 0 AA34 ) ( AA1 * )
A := P' A P = (-----------+-----------) = ( ),
( 0 0 | AA22 0 ) ( 0 AA2 )
( | )
( 0 0 | 0 AA44 )
( 0 BB13 | 0 BB14 )
( | )
( BB31 0 | BB32 0 ) ( BB1 * )
B := P' B P = (-----------+-----------) = ( ).
( 0 0 | 0 BB24 ) ( 0 BB2 )
( | )
( 0 0 | BB42 0 )
Then, further orthogonal transformations that are provided by
MB03FD and MB03DD are used to triangularize the subpencil
aAA1 - bBB1.
Finally, the subpencil aAA2 - bBB2 is triangularized by applying a
special permutation matrix.
See also page 31 in [1] for more details.
</PRE>
<A name="References"><B><FONT SIZE="+1">References</FONT></B></A>
<PRE>
[1] Benner, P., Byers, R., Losse, P., Mehrmann, V. and Xu, H.
Numerical Solution of Real Skew-Hamiltonian/Hamiltonian
Eigenproblems.
Tech. Rep., Technical University Chemnitz, Germany,
Nov. 2007.
[2] Bojanczyk, A., Golub, G. H. and Van Dooren, P.
The periodic Schur decomposition: algorithms and applications.
In F.T. Luk (editor), Advanced Signal Processing Algorithms,
Architectures, and Implementations III, Proc. SPIE Conference,
vol. 1770, pp. 31-42, 1992.
[3] Hench, J. J. and Laub, A. J.
Numerical Solution of the discrete-time periodic Riccati
equation. IEEE Trans. Automat. Control, 39, 1197-1210, 1994.
</PRE>
<A name="Numerical Aspects"><B><FONT SIZE="+1">Numerical Aspects</FONT></B></A>
<PRE> 3
The algorithm is numerically backward stable and needs O(N ) real
floating point operations.
</PRE>
<A name="Comments"><B><FONT SIZE="+1">Further Comments</FONT></B></A>
<PRE>
None
</PRE>
<A name="Example"><B><FONT SIZE="+1">Example</FONT></B></A>
<P>
<B>Program Text</B>
<PRE>
None
</PRE>
<B>Program Data</B>
<PRE>
None
</PRE>
<B>Program Results</B>
<PRE>
None
</PRE>
<HR>
<A HREF=support.html><B>Return to Supporting Routines index</B></A></BODY>
</HTML>
|
