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
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
|
<HTML>
<HEAD><TITLE>TG01LD - SLICOT Library Routine Documentation</TITLE>
</HEAD>
<BODY>
<H2><A Name="TG01LD">TG01LD</A></H2>
<H3>
Finite-infinite decomposition of a descriptor system
</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 orthogonal transformation matrices Q and Z which
reduce the regular pole pencil A-lambda*E of the descriptor system
(A-lambda*E,B,C) to the form (if JOB = 'F')
( Af * ) ( Ef * )
Q'*A*Z = ( ) , Q'*E*Z = ( ) , (1)
( 0 Ai ) ( 0 Ei )
or to the form (if JOB = 'I')
( Ai * ) ( Ei * )
Q'*A*Z = ( ) , Q'*E*Z = ( ) , (2)
( 0 Af ) ( 0 Ef )
where the subpencil Af-lambda*Ef, with Ef nonsingular and upper
triangular, contains the finite eigenvalues, and the subpencil
Ai-lambda*Ei, with Ai nonsingular and upper triangular, contains
the infinite eigenvalues. The subpencil Ai-lambda*Ei is in a
staircase form (see METHOD). If JOBA = 'H', the submatrix Af
is further reduced to an upper Hessenberg form.
</PRE>
<A name="Specification"><B><FONT SIZE="+1">Specification</FONT></B></A>
<PRE>
SUBROUTINE TG01LD( JOB, JOBA, COMPQ, COMPZ, N, M, P, A, LDA,
$ E, LDE, B, LDB, C, LDC, Q, LDQ, Z, LDZ, NF, ND,
$ NIBLCK, IBLCK, TOL, IWORK, DWORK, LDWORK,
$ INFO )
C .. Scalar Arguments ..
CHARACTER COMPQ, COMPZ, JOB, JOBA
INTEGER INFO, LDA, LDB, LDC, LDE, LDQ, LDWORK, LDZ, M,
$ N, ND, NF, NIBLCK, P
DOUBLE PRECISION TOL
C .. Array Arguments ..
INTEGER IBLCK( * ), IWORK(*)
DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ),
$ DWORK( * ), E( LDE, * ), Q( LDQ, * ),
$ Z( LDZ, * )
</PRE>
<A name="Arguments"><B><FONT SIZE="+1">Arguments</FONT></B></A>
<P>
<B>Mode Parameters</B>
<PRE>
JOB CHARACTER*1
= 'F': perform the finite-infinite separation;
= 'I': perform the infinite-finite separation.
JOBA CHARACTER*1
= 'H': reduce Af further to an upper Hessenberg form;
= 'N': keep Af unreduced.
COMPQ CHARACTER*1
= 'N': do not compute Q;
= 'I': Q is initialized to the unit matrix, and the
orthogonal matrix Q is returned;
= 'U': Q must contain an orthogonal matrix Q1 on entry,
and the product Q1*Q is returned.
COMPZ CHARACTER*1
= 'N': do not compute Z;
= 'I': Z is initialized to the unit matrix, and the
orthogonal matrix Z is returned;
= 'U': Z must contain an orthogonal matrix Z1 on entry,
and the product Z1*Z is returned.
</PRE>
<B>Input/Output Parameters</B>
<PRE>
N (input) INTEGER
The number of rows of the matrix B, the number of columns
of the matrix C and the order of the square matrices A
and E. N >= 0.
M (input) INTEGER
The number of columns of the matrix B. M >= 0.
P (input) INTEGER
The number of rows of the matrix C. P >= 0.
A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
On entry, the leading N-by-N part of this array must
contain the N-by-N state matrix A.
On exit, the leading N-by-N part of this array contains
the transformed state matrix Q'*A*Z,
( Af * ) ( Ai * )
Q'*A*Z = ( ) , or Q'*A*Z = ( ) ,
( 0 Ai ) ( 0 Af )
depending on JOB, with Af an NF-by-NF matrix, and Ai an
(N-NF)-by-(N-NF) nonsingular and upper triangular matrix.
If JOBA = 'H', Af is in an upper Hessenberg form.
Otherwise, Af is unreduced.
Ai has a block structure as in (3) or (4), where A0,0 is
ND-by-ND and Ai,i , for i = 1, ..., NIBLCK, is
IBLCK(i)-by-IBLCK(i).
LDA INTEGER
The leading dimension of the array A. LDA >= MAX(1,N).
E (input/output) DOUBLE PRECISION array, dimension (LDE,N)
On entry, the leading N-by-N part of this array must
contain the N-by-N descriptor matrix E.
On exit, the leading N-by-N part of this array contains
the transformed descriptor matrix Q'*E*Z,
( Ef * ) ( Ei * )
Q'*E*Z = ( ) , or Q'*E*Z = ( ) ,
( 0 Ei ) ( 0 Ef )
depending on JOB, with Ef an NF-by-NF nonsingular matrix,
and Ei an (N-NF)-by-(N-NF) nilpotent matrix in an upper
block triangular form, as in (3) or (4).
LDE INTEGER
The leading dimension of the array E. LDE >= MAX(1,N).
B (input/output) DOUBLE PRECISION array, dimension (LDB,K),
where K = M if JOB = 'F', and K = MAX(M,P) if JOB = 'I'.
On entry, the leading N-by-M part of this array must
contain the N-by-M input matrix B.
On exit, the leading N-by-M part of this array contains
the transformed input matrix Q'*B.
LDB INTEGER
The leading dimension of the array B. LDB >= MAX(1,N).
C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
On entry, the leading P-by-N part of this array must
contain the state/output matrix C.
On exit, the leading P-by-N part of this array contains
the transformed matrix C*Z.
LDC INTEGER
The leading dimension of the array C. LDC >= MAX(1,K),
where K = P if JOB = 'F', and K = MAX(M,P) if JOB = 'I'.
Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
If COMPQ = 'N': Q is not referenced.
If COMPQ = 'I': on entry, Q need not be set;
on exit, the leading N-by-N part of this
array contains the orthogonal matrix Q,
where Q' is the product of Householder
transformations applied to A, E, and B on
the left.
If COMPQ = 'U': on entry, the leading N-by-N part of this
array must contain an orthogonal matrix
Q1;
on exit, the leading N-by-N part of this
array contains the orthogonal matrix
Q1*Q.
LDQ INTEGER
The leading dimension of the array Q.
LDQ >= 1, if COMPQ = 'N';
LDQ >= MAX(1,N), if COMPQ = 'I' or 'U'.
Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
If COMPZ = 'N': Z is not referenced.
If COMPZ = 'I': on entry, Z need not be set;
on exit, the leading N-by-N part of this
array contains the orthogonal matrix Z,
which is the product of Householder
transformations applied to A, E, and C on
the right.
If COMPZ = 'U': on entry, the leading N-by-N part of this
array must contain an orthogonal matrix
Z1;
on exit, the leading N-by-N part of this
array contains the orthogonal matrix
Z1*Z.
LDZ INTEGER
The leading dimension of the array Z.
LDZ >= 1, if COMPZ = 'N';
LDZ >= MAX(1,N), if COMPZ = 'I' or 'U'.
NF (output) INTEGER.
The order of the reduced matrices Af and Ef; also, the
number of finite generalized eigenvalues of the pencil
A-lambda*E.
ND (output) INTEGER.
The number of non-dynamic infinite eigenvalues of the
pair (A,E). Note: N-ND is the rank of the matrix E.
NIBLCK (output) INTEGER
If ND > 0, the number of infinite blocks minus one.
If ND = 0, then NIBLCK = 0.
IBLCK (output) INTEGER array, dimension (N)
IBLCK(i) contains the dimension of the i-th block in the
staircase form (3) or (4), with i = 1,2, ..., NIBLCK.
</PRE>
<B>Tolerances</B>
<PRE>
TOL DOUBLE PRECISION
A tolerance used in rank decisions to determine the
effective rank, which is defined as the order of the
largest leading (or trailing) triangular submatrix in the
QR factorization with column pivoting whose estimated
condition number is less than 1/TOL. If the user sets
TOL <= 0, then an implicitly computed, default tolerance,
TOLDEF = N**2*EPS, is used instead, where EPS is the
machine precision (see LAPACK Library routine DLAMCH).
TOL < 1.
</PRE>
<B>Workspace</B>
<PRE>
IWORK INTEGER array, dimension (N)
DWORK DOUBLE PRECISION array, dimension (LDWORK)
On exit, if INFO = 0, DWORK(1) returns the optimal value
of LDWORK.
LDWORK INTEGER
The length of the array DWORK. LDWORK >= 1, and if N > 0,
LDWORK >= N + MAX(3*N,M,P).
If LDWORK = -1, then 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 related to LDWORK is issued by
XERBLA.
</PRE>
<B>Error Indicator</B>
<PRE>
INFO INTEGER
= 0: successful exit;
< 0: if INFO = -i, the i-th argument had an illegal
value;
= 1: the pencil A-lambda*E is not regular.
</PRE>
<A name="Method"><B><FONT SIZE="+1">Method</FONT></B></A>
<PRE>
The subroutine is based on the reduction algorithm of [1].
If JOB = 'F', the matrices Ai and Ei have the form
( A0,0 A0,k ... A0,1 ) ( 0 E0,k ... E0,1 )
Ai = ( 0 Ak,k ... Ak,1 ) , Ei = ( 0 0 ... Ek,1 ) ; (3)
( : : . : ) ( : : . : )
( 0 0 ... A1,1 ) ( 0 0 ... 0 )
if JOB = 'I', the matrices Ai and Ei have the form
( A1,1 ... A1,k A1,0 ) ( 0 ... E1,k E1,0 )
Ai = ( : . : : ) , Ei = ( : . : : ) , (4)
( : ... Ak,k Ak,0 ) ( : ... 0 Ek,0 )
( 0 ... 0 A0,0 ) ( 0 ... 0 0 )
where Ai,i , for i = 0, 1, ..., k, are nonsingular upper
triangular matrices. A0,0 corresponds to the non-dynamic infinite
modes of the system.
</PRE>
<A name="References"><B><FONT SIZE="+1">References</FONT></B></A>
<PRE>
[1] Misra, P., Van Dooren, P., and Varga, A.
Computation of structural invariants of generalized
state-space systems.
Automatica, 30, pp. 1921-1936, 1994.
</PRE>
<A name="Numerical Aspects"><B><FONT SIZE="+1">Numerical Aspects</FONT></B></A>
<PRE>
The algorithm is numerically backward stable and requires
0( N**3 ) floating point operations.
</PRE>
<A name="Comments"><B><FONT SIZE="+1">Further Comments</FONT></B></A>
<PRE>
The number of infinite poles is computed as
NIBLCK
NINFP = Sum IBLCK(i) = N - ND - NF.
i=1
The multiplicities of infinite poles can be computed as follows:
there are IBLCK(k)-IBLCK(k+1) infinite poles of multiplicity
k, for k = 1, ..., NIBLCK, where IBLCK(NIBLCK+1) = 0.
Note that each infinite pole of multiplicity k corresponds to
an infinite eigenvalue of multiplicity k+1.
</PRE>
<A name="Example"><B><FONT SIZE="+1">Example</FONT></B></A>
<P>
<B>Program Text</B>
<PRE>
* TG01LD EXAMPLE PROGRAM TEXT
*
* .. Parameters ..
INTEGER NIN, NOUT
PARAMETER ( NIN = 5, NOUT = 6 )
INTEGER NMAX, MMAX, PMAX
PARAMETER ( NMAX = 20, MMAX = 20, PMAX = 20 )
INTEGER LDA, LDB, LDC, LDE, LDQ, LDZ
PARAMETER ( LDA = NMAX, LDB = NMAX, LDC = PMAX,
$ LDE = NMAX, LDQ = NMAX, LDZ = NMAX )
INTEGER LDWORK
PARAMETER ( LDWORK = NMAX+MAX( 3*NMAX, MMAX, PMAX ) )
* .. Local Scalars ..
CHARACTER*1 COMPQ, COMPZ, JOB, JOBA
INTEGER I, INFO, J, M, N, ND, NF, NIBLCK, P
DOUBLE PRECISION TOL
* .. Local Arrays ..
INTEGER IBLCK(NMAX), IWORK(NMAX)
DOUBLE PRECISION A(LDA,NMAX), B(LDB,MMAX), C(LDC,NMAX),
$ DWORK(LDWORK), E(LDE,NMAX), Q(LDQ,NMAX),
$ Z(LDZ,NMAX)
* .. External Subroutines ..
EXTERNAL TG01LD
* .. Intrinsic Functions ..
INTRINSIC MAX
* .. Executable Statements ..
*
WRITE ( NOUT, FMT = 99999 )
* Skip the heading in the data file and read the data.
READ ( NIN, FMT = '()' )
READ ( NIN, FMT = * ) N, M, P, JOB, JOBA, TOL
COMPQ = 'I'
COMPZ = 'I'
IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99988 ) N
ELSE
READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N )
READ ( NIN, FMT = * ) ( ( E(I,J), J = 1,N ), I = 1,N )
IF ( M.LT.0 .OR. M.GT.MMAX ) THEN
WRITE ( NOUT, FMT = 99987 ) M
ELSE
READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,M ), I = 1,N )
IF ( P.LT.0 .OR. P.GT.PMAX ) THEN
WRITE ( NOUT, FMT = 99986 ) P
ELSE
READ ( NIN, FMT = * ) ( ( C(I,J), J = 1,N ), I = 1,P )
* Find the reduced descriptor system
* (A-lambda E,B,C).
CALL TG01LD( JOB, JOBA, COMPQ, COMPZ, N, M, P, A, LDA,
$ E, LDE, B, LDB, C, LDC, Q, LDQ, Z, LDZ, NF,
$ ND, NIBLCK, IBLCK, TOL, IWORK, DWORK,
$ LDWORK, INFO )
*
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
WRITE ( NOUT, FMT = 99994 ) NF, ND
WRITE ( NOUT, FMT = 99989 ) NIBLCK + 1
IF ( NIBLCK.GT.0 ) THEN
WRITE ( NOUT, FMT = 99985 )
$ ( IBLCK(I), I = 1, NIBLCK )
END IF
WRITE ( NOUT, FMT = 99997 )
DO 10 I = 1, N
WRITE ( NOUT, FMT = 99995 ) ( A(I,J), J = 1,N )
10 CONTINUE
WRITE ( NOUT, FMT = 99996 )
DO 20 I = 1, N
WRITE ( NOUT, FMT = 99995 ) ( E(I,J), J = 1,N )
20 CONTINUE
WRITE ( NOUT, FMT = 99993 )
DO 30 I = 1, N
WRITE ( NOUT, FMT = 99995 ) ( B(I,J), J = 1,M )
30 CONTINUE
WRITE ( NOUT, FMT = 99992 )
DO 40 I = 1, P
WRITE ( NOUT, FMT = 99995 ) ( C(I,J), J = 1,N )
40 CONTINUE
WRITE ( NOUT, FMT = 99991 )
DO 50 I = 1, N
WRITE ( NOUT, FMT = 99995 ) ( Q(I,J), J = 1,N )
50 CONTINUE
WRITE ( NOUT, FMT = 99990 )
DO 60 I = 1, N
WRITE ( NOUT, FMT = 99995 ) ( Z(I,J), J = 1,N )
60 CONTINUE
END IF
END IF
END IF
END IF
STOP
*
99999 FORMAT (' TG01LD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from TG01LD = ',I2)
99997 FORMAT (/' The reduced state dynamics matrix Q''*A*Z is ')
99996 FORMAT (/' The reduced descriptor matrix Q''*E*Z is ')
99995 FORMAT (20(1X,F8.4))
99994 FORMAT (' Order of reduced system =', I5/
$ ' Number of non-dynamic infinite eigenvalues =', I5)
99993 FORMAT (/' The reduced input/state matrix Q''*B is ')
99992 FORMAT (/' The reduced state/output matrix C*Z is ')
99991 FORMAT (/' The left transformation matrix Q is ')
99990 FORMAT (/' The right transformation matrix Z is ')
99989 FORMAT ( ' Number of infinite blocks = ',I5)
99988 FORMAT (/' N is out of range.',/' N = ',I5)
99987 FORMAT (/' M is out of range.',/' M = ',I5)
99986 FORMAT (/' P is out of range.',/' P = ',I5)
99985 FORMAT ( ' Dimension of the blocks'/20I5)
END
</PRE>
<B>Program Data</B>
<PRE>
TG01LD EXAMPLE PROGRAM DATA
4 2 2 F N 0.0
-1 0 0 3
0 0 1 2
1 1 0 4
0 0 0 0
1 2 0 0
0 1 0 1
3 9 6 3
0 0 2 0
1 0
0 0
0 1
1 1
-1 0 1 0
0 1 -1 1
</PRE>
<B>Program Results</B>
<PRE>
TG01LD EXAMPLE PROGRAM RESULTS
Order of reduced system = 3
Number of non-dynamic infinite eigenvalues = 1
Number of infinite blocks = 1
The reduced state dynamics matrix Q'*A*Z is
2.4497 -1.3995 0.2397 -4.0023
-0.0680 -0.0030 0.1739 -1.6225
0.3707 0.0161 -0.9482 0.1049
0.0000 0.0000 0.0000 2.2913
The reduced descriptor matrix Q'*E*Z is
9.9139 4.7725 -3.4725 -2.3836
0.0000 -1.2024 2.0137 0.7926
0.0000 0.0000 0.2929 -0.9914
0.0000 0.0000 0.0000 0.0000
The reduced input/state matrix Q'*B is
-0.2157 -0.9705
0.3015 0.9516
0.7595 0.0991
1.1339 0.3780
The reduced state/output matrix C*Z is
0.5345 -1.1134 0.3758 0.5774
-1.0690 0.2784 -1.2026 0.5774
The left transformation matrix Q is
-0.2157 -0.5088 0.6109 0.5669
-0.1078 -0.2544 -0.7760 0.5669
-0.9705 0.1413 -0.0495 -0.1890
0.0000 0.8102 0.1486 0.5669
The right transformation matrix Z is
-0.5345 0.6263 0.4617 -0.3299
-0.8018 -0.5219 -0.2792 -0.0825
0.0000 -0.4871 0.8375 0.2474
-0.2673 0.3132 -0.0859 0.9073
</PRE>
<HR>
<p>
<A HREF=..\libindex.html><B>Return to index</B></A></BODY>
</HTML>
|
